[Scipy-svn] r2607 - in trunk/Lib/linalg: . tests

scipy-svn at scipy.org scipy-svn at scipy.org
Wed Jan 24 08:35:05 CST 2007


Author: edschofield
Date: 2007-01-24 08:35:01 -0600 (Wed, 24 Jan 2007)
New Revision: 2607

Added:
   trunk/Lib/linalg/tests/test_iterative.py
Modified:
   trunk/Lib/linalg/iterative.py
Log:
Added callbacks to linalg.iterative solvers


Modified: trunk/Lib/linalg/iterative.py
===================================================================
--- trunk/Lib/linalg/iterative.py	2007-01-24 13:53:11 UTC (rev 2606)
+++ trunk/Lib/linalg/iterative.py	2007-01-24 14:35:01 UTC (rev 2607)
@@ -108,7 +108,7 @@
 class get_rpsolveq(get_psolveq):
     methname = 'rpsolve'
 
-def bicg(A,b,x0=None,tol=1e-5,maxiter=None,xtype=None):
+def bicg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
     """Use BIConjugate Gradient iteration to solve A x = b
 
     Inputs:
@@ -133,11 +133,15 @@
     x0  -- (0) default starting guess
     tol -- (1e-5) relative tolerance to achieve
     maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be determined
-                 from A.dtype.char and b.  If A does not have a typecode method then it will
-                 compute A.matvec(x0) to get a typecode.   To save the extra computation when
-                 A does not have a typecode attribute use xtype=0 for the same type as b or
-                 use xtype='f','d','F',or 'D'
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
     """
     b = sb.asarray(b)+0.0
     n = len(b)
@@ -180,9 +184,12 @@
     ftflag = True
     bnrm2 = -1.0
     iter_ = maxiter
-    while 1:
+    while True:
+        olditer = iter_
         x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
            revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
         slice1 = slice(ndx1-1, ndx1-1+n)
         slice2 = slice(ndx2-1, ndx2-1+n)
         if (ijob == -1):
@@ -219,7 +226,8 @@
 
     return x, info
 
-def bicgstab(A,b,x0=None,tol=1e-5,maxiter=None,xtype=None):
+
+def bicgstab(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
     """Use BIConjugate Gradient STABilized iteration to solve A x = b
 
     Inputs:
@@ -243,11 +251,15 @@
     x0  -- (0) default starting guess
     tol -- (1e-5) relative tolerance to achieve
     maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be determined
-                 from A.dtype.char and b.  If A does not have a typecode method then it will
-                 compute A.matvec(x0) to get a typecode.   To save the extra computation when
-                 A does not have a typecode attribute use xtype=0 for the same type as b or
-                 use xtype='f','d','F',or 'D'
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
     """
     b = sb.asarray(b)+0.0
     n = len(b)
@@ -290,9 +302,12 @@
     ftflag = True
     bnrm2 = -1.0
     iter_ = maxiter
-    while 1:
+    while True:
+        olditer = iter_
         x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
            revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
         slice1 = slice(ndx1-1, ndx1-1+n)
         slice2 = slice(ndx2-1, ndx2-1+n)
         if (ijob == -1):
@@ -320,7 +335,8 @@
 
     return x, info
 
-def cg(A,b,x0=None,tol=1e-5,maxiter=None,xtype=None):
+
+def cg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
     """Use Conjugate Gradient iteration to solve A x = b (A^H = A)
 
     Inputs:
@@ -345,11 +361,15 @@
     x0  -- (0) default starting guess
     tol -- (1e-5) relative tolerance to achieve
     maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be determined
-                 from A.dtype.char and b.  If A does not have a typecode method then it will
-                 compute A.matvec(x0) to get a typecode.   To save the extra computation when
-                 A does not have a typecode attribute use xtype=0 for the same type as b or
-                 use xtype='f','d','F',or 'D'
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
     """
     b = sb.asarray(b)+0.0
     n = len(b)
@@ -392,9 +412,12 @@
     ftflag = True
     bnrm2 = -1.0
     iter_ = maxiter
-    while 1:
+    while True:
+        olditer = iter_
         x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
            revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
         slice1 = slice(ndx1-1, ndx1-1+n)
         slice2 = slice(ndx2-1, ndx2-1+n)
         if (ijob == -1):
@@ -423,7 +446,7 @@
     return x, info
 
 
-def cgs(A,b,x0=None,tol=1e-5,maxiter=None,xtype=None):
+def cgs(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
     """Use Conjugate Gradient Squared iteration to solve A x = b
 
     Inputs:
@@ -448,11 +471,15 @@
     x0  -- (0) default starting guess
     tol -- (1e-5) relative tolerance to achieve
     maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be determined
-                 from A.dtype.char and b.  If A does not have a typecode method then it will
-                 compute A.matvec(x0) to get a typecode.   To save the extra computation when
-                 A does not have a typecode attribute use xtype=0 for the same type as b or
-                 use xtype='f','d','F',or 'D'
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
     """
     b = sb.asarray(b) + 0.0
     n = len(b)
@@ -495,9 +522,12 @@
     ftflag = True
     bnrm2 = -1.0
     iter_ = maxiter
-    while 1:
+    while True:
+        olditer = iter_
         x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
            revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
         slice1 = slice(ndx1-1, ndx1-1+n)
         slice2 = slice(ndx2-1, ndx2-1+n)
         if (ijob == -1):
@@ -525,7 +555,8 @@
 
     return x, info
 
-def gmres(A,b,restrt=None,x0=None,tol=1e-5,maxiter=None,xtype=None):
+
+def gmres(A, b, x0=None, tol=1e-5, restrt=None, maxiter=None, xtype=None, callback=None):
     """Use Generalized Minimal RESidual iteration to solve A x = b
 
     Inputs:
@@ -552,11 +583,15 @@
     x0  -- (0) default starting guess
     tol -- (1e-5) relative tolerance to achieve
     maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be determined
-                 from A.dtype.char and b.  If A does not have a dtype attribute then it will
-                 compute A.matvec(x0) to get a data type.   To save the extra computation when
-                 A does not have a typecode attribute use xtype=0 for the same type as b or
-                 use xtype='f','d','F',or 'D'
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
     """
     b = sb.asarray(b)+0.0
     n = len(b)
@@ -600,9 +635,12 @@
     ftflag = True
     bnrm2 = -1.0
     iter_ = maxiter
-    while 1:
+    while True:
+        olditer = iter_
         x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
            revcom(b, x, restrt, work, work2, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
         slice1 = slice(ndx1-1, ndx1-1+n)
         slice2 = slice(ndx2-1, ndx2-1+n)
         if (ijob == -1):
@@ -631,8 +669,8 @@
     return x, info
 
 
-def qmr(A,b,x0=None,tol=1e-5,maxiter=None,xtype=None):
-    """Use Quasi-Minimal Residucal iteration to solve A x = b
+def qmr(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None):
+    """Use Quasi-Minimal Residual iteration to solve A x = b
 
     Inputs:
 
@@ -657,11 +695,15 @@
     x0  -- (0) default starting guess
     tol -- (1e-5) relative tolerance to achieve
     maxiter -- (10*n) maximum number of iterations
-    xtype  --  The type of the result.  If None, then it will be determined
-                 from A.dtype.char and b.  If A does not have a typecode method then it will
-                 compute A.matvec(x0) to get a typecode.   To save the extra computation when
-                 A does not have a typecode attribute use xtype=0 for the same type as b or
-                 use xtype='f','d','F',or 'D'
+    xtype  --  The type of the result.  If None, then it will be
+               determined from A.dtype.char and b.  If A does not have a
+               typecode method then it will compute A.matvec(x0) to get a
+               typecode.   To save the extra computation when A does not
+               have a typecode attribute use xtype=0 for the same type as
+               b or use xtype='f','d','F',or 'D'
+    callback -- an optional user-supplied function to call after each
+                iteration.  It is called as callback(xk), where xk is the
+                current parameter vector.
     """
     b = sb.asarray(b)+0.0
     n = len(b)
@@ -705,9 +747,12 @@
     ftflag = True
     bnrm2 = -1.0
     iter_ = maxiter
-    while 1:
+    while True:
+        olditer = iter_
         x, iter_, resid, info, ndx1, ndx2, sclr1, sclr2, ijob = \
            revcom(b, x, work, iter_, resid, info, ndx1, ndx2, ijob)
+        if callback is not None and iter_ > olditer:
+            callback(x)
         slice1 = slice(ndx1-1, ndx1-1+n)
         slice2 = slice(ndx2-1, ndx2-1+n)
         if (ijob == -1):

Added: trunk/Lib/linalg/tests/test_iterative.py
===================================================================
--- trunk/Lib/linalg/tests/test_iterative.py	2007-01-24 13:53:11 UTC (rev 2606)
+++ trunk/Lib/linalg/tests/test_iterative.py	2007-01-24 14:35:01 UTC (rev 2607)
@@ -0,0 +1,73 @@
+#!/usr/bin/env python
+#
+# Created by: Ed Schofield, Jan 2007
+#
+""" Test functions for the linalg.iterative module
+"""
+__usage__ = """
+Build linalg:
+  python setup_linalg.py build
+Run tests if scipy is installed:
+  python -c 'import scipy;scipy.linalg.test(<level>)'
+Run tests if linalg is not installed:
+  python tests/test_iterative.py [<level>]
+"""
+
+from numpy import zeros, dot, diag, ones
+from numpy.testing import *
+from numpy.random import rand
+#from numpy import arange, add, array, dot, zeros, identity, conjugate, transpose
+
+import sys
+set_package_path()
+from linalg import iterative, norm, cg, cgs, bicg, bicgstab, gmres, qmr
+restore_path()
+
+
+def callback(x):
+    global A, b
+    res = b-dot(A,x)
+    print "||A.x - b|| = " + str(norm(dot(A,x)-b))
+
+class test_iterative_solvers(NumpyTestCase):
+    def __init__(self, *args, **kwds):
+        NumpyTestCase.__init__(self, *args, **kwds)
+        self.setUp()
+    def setUp (self):
+        global A, b
+        n = 10
+        self.tol = 1e-5
+        self.x0 = zeros(n, float)
+        A = rand(n, n)+diag(4*ones(n))
+        self.A = 0.5 * (A+A.T)
+        A = self.A
+        self.b = rand(n)
+        b = self.b
+
+    def check_cg(self):
+        x, info = cg(self.A, self.b, self.x0, callback=callback)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def check_bicg(self):
+        x, info = bicg(self.A, self.b, self.x0, callback=callback)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def check_cgs(self):
+        x, info = cgs(self.A, self.b, self.x0, callback=callback)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def check_bicgstab(self):
+        x, info = bicgstab(self.A, self.b, self.x0, callback=callback)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def check_gmres(self):
+        x, info = gmres(self.A, self.b, self.x0, callback=callback)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+    def check_qmr(self):
+        x, info = qmr(self.A, self.b, self.x0, callback=callback)
+        assert norm(dot(self.A, x) - self.b) < 5*self.tol
+
+if __name__ == "__main__":
+    NumpyTest().run()
+



More information about the Scipy-svn mailing list