[Scipy-svn] r5034 - in trunk/scipy/sparse/linalg/isolve: . tests

scipy-svn@scip... scipy-svn@scip...
Sun Nov 9 21:06:01 CST 2008


Author: wnbell
Date: 2008-11-09 21:05:57 -0600 (Sun, 09 Nov 2008)
New Revision: 5034

Modified:
   trunk/scipy/sparse/linalg/isolve/iterative.py
   trunk/scipy/sparse/linalg/isolve/minres.py
   trunk/scipy/sparse/linalg/isolve/tests/test_iterative.py
Log:
test 'info' output of iterative methods
make GMRES use default restart of 10 (as opposed to n)
resolves ticket #666 (the mark of the beast!)


Modified: trunk/scipy/sparse/linalg/isolve/iterative.py
===================================================================
--- trunk/scipy/sparse/linalg/isolve/iterative.py	2008-11-10 01:57:55 UTC (rev 5033)
+++ trunk/scipy/sparse/linalg/isolve/iterative.py	2008-11-10 03:05:57 UTC (rev 5034)
@@ -1,14 +1,5 @@
-## Automatically adapted for scipy Oct 18, 2005 by
+"""Iterative methods for solving linear systems"""
 
-
-# Iterative methods using reverse-communication raw material
-#   These methods solve
-#   Ax = b  for x
-
-#   where A must have A.matvec(x,*args) defined
-#    or be a numeric array
-
-
 __all__ = ['bicg','bicgstab','cg','cgs','gmres','qmr']
 
 import _iterative
@@ -106,6 +97,10 @@
                 ftflag = False
             bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
         ijob = 2
+    
+    if info > 0 and iter_ == maxiter and resid > tol:
+        #info isn't set appropriately otherwise
+        info = iter_
 
     return postprocess(x), info
 
@@ -197,6 +192,10 @@
                 ftflag = False
             bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
         ijob = 2
+    
+    if info > 0 and iter_ == maxiter and resid > tol:
+        #info isn't set appropriately otherwise
+        info = iter_
 
     return postprocess(x), info
 
@@ -284,6 +283,11 @@
             bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
         ijob = 2
 
+
+    if info > 0 and iter_ == maxiter and resid > tol:
+        #info isn't set appropriately otherwise
+        info = iter_
+
     return postprocess(x), info
 
 
@@ -369,11 +373,15 @@
                 ftflag = False
             bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
         ijob = 2
+    
+    if info > 0 and iter_ == maxiter and resid > tol:
+        #info isn't set appropriately otherwise
+        info = iter_
 
     return postprocess(x), info
 
 
-def gmres(A, b, x0=None, tol=1e-5, restrt=None, maxiter=None, xtype=None, M=None, callback=None):
+def gmres(A, b, x0=None, tol=1e-5, restrt=20, maxiter=None, xtype=None, M=None, callback=None):
     """Use Generalized Minimal RESidual iteration to solve A x = b
 
     Inputs:
@@ -397,7 +405,7 @@
 
     x0  -- (0) default starting guess.
     tol -- (1e-5) relative tolerance to achieve
-    restrt -- (n) When to restart (change this to get faster performance -- but
+    restrt -- (10) When to restart (change this to get faster performance -- but
                    may not converge).
     maxiter -- (10*n) maximum number of iterations
     xtype  --  The type of the result.  If None, then it will be
@@ -416,14 +424,14 @@
     if maxiter is None:
         maxiter = n*10
 
+    restrt = min(restrt, n)        
+
     matvec = A.matvec
     psolve = M.matvec
     ltr = _type_conv[x.dtype.char]
     revcom   = getattr(_iterative, ltr + 'gmresrevcom')
     stoptest = getattr(_iterative, ltr + 'stoptest2')
 
-    if restrt is None:
-        restrt = n
     resid = tol
     ndx1 = 1
     ndx2 = -1
@@ -480,6 +488,10 @@
 
         if iter_num > maxiter:
             break
+    
+    if info >= 0 and resid > tol:
+        #info isn't set appropriately otherwise
+        info = maxiter
 
     return postprocess(x), info
 
@@ -593,6 +605,10 @@
                 ftflag = False
             bnrm2, resid, info = stoptest(work[slice1], b, bnrm2, tol, info)
         ijob = 2
+    
+    if info > 0 and iter_ == maxiter and resid > tol:
+        #info isn't set appropriately otherwise
+        info = iter_
 
     return postprocess(x), info
 

Modified: trunk/scipy/sparse/linalg/isolve/minres.py
===================================================================
--- trunk/scipy/sparse/linalg/isolve/minres.py	2008-11-10 01:57:55 UTC (rev 5033)
+++ trunk/scipy/sparse/linalg/isolve/minres.py	2008-11-10 03:05:57 UTC (rev 5034)
@@ -67,7 +67,7 @@
     istop = 0;   itn   = 0;   Anorm = 0;    Acond = 0;
     rnorm = 0;   ynorm = 0;
 
-    xtype = A.dtype #TODO update
+    xtype = x.dtype
 
     eps = finfo(xtype).eps
 
@@ -273,9 +273,14 @@
         print last + ' Arnorm  =  %12.4e'                       %  (Arnorm,)
         print last + msg[istop+1]
 
-    return (postprocess(x),0)
+    if istop == 6:
+        info = maxiter
+    else:
+        info = 0
 
+    return (postprocess(x),info)
 
+
 if __name__ == '__main__':
     from scipy import ones, arange
     from scipy.linalg import norm

Modified: trunk/scipy/sparse/linalg/isolve/tests/test_iterative.py
===================================================================
--- trunk/scipy/sparse/linalg/isolve/tests/test_iterative.py	2008-11-10 01:57:55 UTC (rev 5033)
+++ trunk/scipy/sparse/linalg/isolve/tests/test_iterative.py	2008-11-10 03:05:57 UTC (rev 5034)
@@ -8,6 +8,7 @@
 from scipy.linalg import norm
 from scipy.sparse import spdiags, csr_matrix
 
+from scipy.sparse.linalg.interface import LinearOperator
 from scipy.sparse.linalg.isolve import cg, cgs, bicg, bicgstab, gmres, qmr, minres
 
 #def callback(x):
@@ -66,6 +67,7 @@
         """test whether maxiter is respected"""
 
         A = Poisson1D
+        tol = 1e-12
 
         for solver,req_sym,req_pos in self.solvers:
             b  = arange(A.shape[0], dtype=float)
@@ -75,11 +77,11 @@
             def callback(x):
                 residuals.append( norm(b - A*x) )
 
-            x, info = solver(A, b, x0=x0, tol=1e-8, maxiter=3, callback=callback)
+            x, info = solver(A, b, x0=x0, tol=tol, maxiter=3, callback=callback)
 
             assert_equal(len(residuals), 3)
+            assert_equal(info, 3)
 
-
     def test_convergence(self):
         """test whether all methods converge"""
 
@@ -101,33 +103,44 @@
                 assert( norm(b - A*x) < tol*norm(b) )
 
     def test_precond(self):
-        """test whether all methods accept a preconditioner"""
+        """test whether all methods accept a trivial preconditioner"""
 
         tol = 1e-8
+        
+        def identity(b,which=None):
+            """trivial preconditioner"""
+            return b
 
         for solver,req_sym,req_pos in self.solvers:
+
             for A,sym,pos in self.cases:
                 if req_sym and not sym: continue
                 if req_pos and not pos: continue
 
                 M,N = A.shape
-                D = spdiags( [abs(1.0/A.diagonal())], [0], M, N)
-                def precond(b,which=None):
-                    return D*b
+                D = spdiags( [1.0/A.diagonal()], [0], M, N)
 
-                A = A.copy()
-                A.psolve  = precond
-                A.rpsolve = precond
-
                 b  = arange(A.shape[0], dtype=float)
                 x0 = 0*b
 
-                x, info = solver(A, b, x0=x0, tol=tol)
+                precond = LinearOperator(A.shape, identity, rmatvec=identity)
 
+                if solver == qmr:
+                    x, info = solver(A, b, M1=precond, M2=precond, x0=x0, tol=tol)
+                else:
+                    x, info = solver(A, b, M=precond, x0=x0, tol=tol)
                 assert_equal(info,0)
                 assert( norm(b - A*x) < tol*norm(b) )
+                
+                A = A.copy()
+                A.psolve  = identity 
+                A.rpsolve = identity
 
+                x, info = solver(A, b, x0=x0, tol=tol)
+                assert_equal(info,0)
+                assert( norm(b - A*x) < tol*norm(b) )
 
+
 class TestQMR(TestCase):
     def test_leftright_precond(self):
         """Check that QMR works with left and right preconditioners"""



More information about the Scipy-svn mailing list