[Scipy-svn] r3888 - in trunk/scipy/splinalg/isolve: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Feb 1 17:29:41 CST 2008


Author: wnbell
Date: 2008-02-01 17:29:37 -0600 (Fri, 01 Feb 2008)
New Revision: 3888

Added:
   trunk/scipy/splinalg/isolve/utils.py
Modified:
   trunk/scipy/splinalg/isolve/__init__.py
   trunk/scipy/splinalg/isolve/minres.py
   trunk/scipy/splinalg/isolve/tests/test_iterative.py
Log:
updated MINRES code
abstracted iterative solver setup code


Modified: trunk/scipy/splinalg/isolve/__init__.py
===================================================================
--- trunk/scipy/splinalg/isolve/__init__.py	2008-02-01 19:49:30 UTC (rev 3887)
+++ trunk/scipy/splinalg/isolve/__init__.py	2008-02-01 23:29:37 UTC (rev 3888)
@@ -2,6 +2,7 @@
 
 #from info import __doc__
 from iterative import *
+from minres import minres
 
 __all__ = filter(lambda s:not s.startswith('_'),dir())
 from scipy.testing.pkgtester import Tester

Modified: trunk/scipy/splinalg/isolve/minres.py
===================================================================
--- trunk/scipy/splinalg/isolve/minres.py	2008-02-01 19:49:30 UTC (rev 3887)
+++ trunk/scipy/splinalg/isolve/minres.py	2008-02-01 23:29:37 UTC (rev 3888)
@@ -1,11 +1,10 @@
-from numpy import sqrt, inner, finfo, asarray, zeros
+from numpy import ndarray, matrix, sqrt, inner, finfo, asarray, zeros
 from numpy.linalg import norm
 
-def psolve(x): return x
-def check_sizes(A,x,b): pass
+from utils import make_system
 
 def minres(A, b, x0=None, shift=0.0, tol=1e-5, maxiter=None, xtype=None,
-        precond=None, callback=None, show=False, check=False):
+           M=None, callback=None, show=False, check=True):
     """Use the Minimum Residual Method (MINRES) to solve Ax=b 
     
     MINRES minimizes norm(A*x - b) for the symmetric matrix A.  Unlike
@@ -30,23 +29,19 @@
             http://www.stanford.edu/group/SOL/software/minres/matlab/
 
     """ 
+    A,M,x,b,postprocess = make_system(A,M,x0,b,xtype)
 
-    show  = True  #TODO remove
-    check = True  #TODO remove
+    matvec = A.matvec
+    psolve = M.matvec
 
     first = 'Enter minres.   '
     last  = 'Exit  minres.   '
 
-    assert(A.shape[0] == A.shape[1])
-    assert(A.shape[1] == len(b))
-
-    b = asarray(b).ravel()
     n = A.shape[0]
 
     if maxiter is None:
         maxiter = 5 * n
 
-    matvec = A.matvec
 
     msg   =[' beta2 = 0.  If M = I, b and x are eigenvectors    ',   # -1
             ' beta1 = 0.  The exact solution is  x = 0          ',   #  0
@@ -56,9 +51,9 @@
             ' x has converged to an eigenvector                 ',   #  4
             ' acond has exceeded 0.1/eps                        ',   #  5
             ' The iteration limit was reached                   ',   #  6
-            ' Aname  does not define a symmetric matrix         ',   #  7
-            ' Mname  does not define a symmetric matrix         ',   #  8
-            ' Mname  does not define a pos-def preconditioner   ']   #  9
+            ' A  does not define a symmetric matrix             ',   #  7
+            ' M  does not define a symmetric matrix             ',   #  8
+            ' M  does not define a pos-def preconditioner       ']   #  9
 
      
     if show:
@@ -90,7 +85,7 @@
     if beta1 < 0:
         raise ValueError('indefinite preconditioner')
     elif beta1 == 0:
-        return x
+        return (postprocess(x), 0)
     
     beta1 = sqrt( beta1 )
 
@@ -262,7 +257,7 @@
         if callback is not None:
             callback(x)
 
-        if istop > 0: break
+        if istop != 0: break #TODO check this
         
 
     if show:
@@ -273,7 +268,7 @@
         print last + ' Arnorm  =  %12.4e'                       %  (Arnorm,)
         print last + msg[istop+1]
 
-    return x
+    return (postprocess(x),0)
 
 
 if __name__ == '__main__':
@@ -283,7 +278,7 @@
     from scipy.splinalg import cg
     #from scipy.sandbox.multigrid import *
 
-    n = 100
+    n = 10
 
     residuals = []
 
@@ -292,7 +287,9 @@
 
     #A = poisson((10,),format='csr')
     A = spdiags( [arange(1,n+1,dtype=float)], [0], n, n, format='csr')
-    b = ones( A.shape[0] )
+    M = spdiags( [1.0/arange(1,n+1,dtype=float)], [0], n, n, format='csr')
+    A.psolve = M.matvec
+    b = 0*ones( A.shape[0] )
     x = minres(A,b,tol=1e-12,maxiter=None,callback=cb)
     #x = cg(A,b,x0=b,tol=1e-12,maxiter=None,callback=cb)[0]
 

Modified: trunk/scipy/splinalg/isolve/tests/test_iterative.py
===================================================================
--- trunk/scipy/splinalg/isolve/tests/test_iterative.py	2008-02-01 19:49:30 UTC (rev 3887)
+++ trunk/scipy/splinalg/isolve/tests/test_iterative.py	2008-02-01 23:29:37 UTC (rev 3888)
@@ -9,7 +9,7 @@
 from scipy.linalg import norm
 from scipy.sparse import spdiags
 
-from scipy.splinalg.isolve import cg, cgs, bicg, bicgstab, gmres, qmr
+from scipy.splinalg.isolve import cg, cgs, bicg, bicgstab, gmres, qmr, minres
 
 #def callback(x):
 #    global A, b
@@ -36,7 +36,7 @@
         self.solvers.append( (bicgstab, False, False) )
         self.solvers.append( (gmres,    False, False) )
         self.solvers.append( (qmr,      False, False) )
-        #self.solvers.append( (minres,   True,  False) )
+        self.solvers.append( (minres,   True,  False) )
         
         # list of tuples (A, symmetric, positive_definite )
         self.cases = []
@@ -91,7 +91,7 @@
                 if req_pos and not pos: continue
 
                 M,N = A.shape
-                D = spdiags( [1.0/A.diagonal()], [0], M, N)
+                D = spdiags( [abs(1.0/A.diagonal())], [0], M, N)
                 def precond(b,which=None):
                     return D*b
 

Added: trunk/scipy/splinalg/isolve/utils.py
===================================================================
--- trunk/scipy/splinalg/isolve/utils.py	2008-02-01 19:49:30 UTC (rev 3887)
+++ trunk/scipy/splinalg/isolve/utils.py	2008-02-01 23:29:37 UTC (rev 3888)
@@ -0,0 +1,79 @@
+from numpy import asanyarray, asmatrix, array, matrix, zeros
+
+from scipy.splinalg.interface import aslinearoperator, LinearOperator
+
+_coerce_rules = {('f','f'):'f', ('f','d'):'d', ('f','F'):'F',
+                 ('f','D'):'D', ('d','f'):'d', ('d','d'):'d',
+                 ('d','F'):'D', ('d','D'):'D', ('F','f'):'F',
+                 ('F','d'):'D', ('F','F'):'F', ('F','D'):'D',
+                 ('D','f'):'D', ('D','d'):'D', ('D','F'):'D',
+                 ('D','D'):'D'}
+
+def coerce(x,y):
+    if x not in 'fdFD':
+        x = 'd'
+    if y not in 'fdFD':
+        y = 'd'
+    return _coerce_rules[x,y]
+
+def id(x):
+    return x
+
+def make_system(A, M, x0, b, xtype=None):
+    A_ = A
+    A = aslinearoperator(A)
+
+    if A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix (shape=%s)' % shape)
+
+    N = A.shape[0]
+    
+    b = asanyarray(b)
+
+    if not (b.shape == (N,1) or b.shape == (N,)):
+        raise ValueError('A and b have incompatible dimensions')
+
+    def postprocess(x):
+        if isinstance(b,matrix):
+            x = asmatrix(x)
+        return x.reshape(b.shape)
+
+
+    if xtype is None:
+        if hasattr(A,'dtype'):
+            xtype = A.dtype.char
+        else:
+            xtype = A.matvec(b).dtype.char
+        xtype = coerce(xtype, b.dtype.char)
+    elif xtype == 0:
+        xtype = b.dtype.char
+    else:
+        if xtype not in 'fdFD':
+            raise ValueError, "xtype must be 'f', 'd', 'F', or 'D'"
+
+    if x0 is None:
+        x = zeros(N, dtype=xtype)
+    else:
+        x = array(x0, dtype=xtype)
+        if not (x.shape == (N,1) or x.shape == (N,)):
+            raise ValueError('A and x have incompatible dimensions')
+        x = x.ravel()
+
+    # process preconditioner
+    if M is None:
+        if hasattr(A_,'psolve'):
+            psolve = A_.psolve
+        else:
+            psolve = id
+        if hasattr(A_,'rpsolve'):
+            rpsolve = A_.rpsolve
+        else:
+            rpsolve = id
+        M = LinearOperator(A.shape, matvec=psolve, rmatvec=rpsolve, dtype=A.dtype)
+    else:
+        if A.shape != M.shape:
+            raise ValueError('matrix and preconditioner have different shapes')
+
+    return A, M, x, b, postprocess
+
+



More information about the Scipy-svn mailing list