[Scipy-svn] r2607 - in trunk/Lib/linalg: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Jan 24 09:35:05 EST 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