[Scipy-svn] r6270 - trunk/scipy/sparse/linalg/eigen/arpack

scipy-svn@scip... scipy-svn@scip...
Fri Mar 26 00:35:07 CDT 2010


Author: cdavid
Date: 2010-03-26 00:35:07 -0500 (Fri, 26 Mar 2010)
New Revision: 6270

Modified:
   trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
Log:
REF: abstrace solver call for ARPACK.

Modified: trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-03-26 05:34:57 UTC (rev 6269)
+++ trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-03-26 05:35:07 UTC (rev 6270)
@@ -52,7 +52,7 @@
 _ndigits = {'f':5, 'd':12, 'F':5, 'D':12}
 
 class _ArpackParams(object):
-    def __init__(self, n, k, tp, mode="symmetric", sigma=None,
+    def __init__(self, n, k, tp, matvec, mode="symmetric", sigma=None,
                  ncv=None, v0=None, maxiter=None, which="LM", tol=0):
         if k <= 0:
             raise ValueError("k must be positive, k=%d" % k)
@@ -95,7 +95,8 @@
 
             self.workd = np.zeros(3 * n, tp)
             self.workl = np.zeros(3 * ncv * ncv + 6 * ncv, tp)
-            self.solver = _arpack.__dict__[ltr + 'naupd']
+            self._arpack_solver = _arpack.__dict__[ltr + 'naupd']
+            self.iterate = self._unsymmetric_solver
             self.extract = _arpack.__dict__[ltr + 'neupd']
 
             if tp in 'FD':
@@ -108,7 +109,8 @@
 
             self.workd = np.zeros(3 * n, tp)
             self.workl = np.zeros(ncv * (ncv + 8), tp)
-            self.solver = _arpack.__dict__[ltr + 'saupd']
+            self._arpack_solver = _arpack.__dict__[ltr + 'saupd']
+            self.iterate = self._symmetric_solver
             self.extract = _arpack.__dict__[ltr + 'seupd']
 
             self.ipntr = np.zeros(11, "int")
@@ -127,6 +129,7 @@
 
         self.n = n
         self.mode = mode
+        self.matvec = matvec
         self.tol = tol
         self.k = k
         self.maxiter = maxiter
@@ -136,6 +139,62 @@
         self.info = info
         self.bmat = 'I'
 
+        self.converged = False
+        self.ido = 0
+
+    def _unsymmetric_solver(self):
+        if self.tp in 'fd':
+            self.ido, self.resid, self.v, self.iparam, self.ipntr, self.info = \
+                self._arpack_solver(self.ido, self.bmat, self.which, self.k, self.tol,
+                        self.resid, self.v, self.iparam, self.ipntr,
+                        self.workd, self.workl, self.info)
+        else:
+            self.ido, self.resid, self.v, self.iparam, self.ipntr, self.info =\
+                self._arpack_solver(self.ido, self.bmat, self.which, self.k, self.tol,
+                        self.resid, self.v, self.iparam, self.ipntr,
+                        self.workd, self.workl, self.rwork, self.info)
+
+        xslice = slice(self.ipntr[0]-1, self.ipntr[0]-1+self.n)
+        yslice = slice(self.ipntr[1]-1, self.ipntr[1]-1+self.n)
+        if self.ido == -1:
+            # initialization
+            self.workd[yslice] = self.matvec(self.workd[xslice])
+        elif self.ido == 1:
+            # compute y=Ax
+            self.workd[yslice] = self.matvec(self.workd[xslice])
+        else:
+            self.converged = True
+
+            if self.info < -1 :
+                raise RuntimeError("Error info=%d in arpack" % self.info)
+            elif self.info == -1:
+                warnings.warn("Maximum number of iterations taken: %s" % self.iparam[2])
+
+    def _symmetric_solver(self):
+        self.ido, self.resid, self.v, self.iparam, self.ipntr, self.info = \
+            self._arpack_solver(self.ido, self.bmat, self.which, self.k, self.tol,
+                    self.resid, self.v, self.iparam, self.ipntr,
+                    self.workd, self.workl, self.info)
+
+        xslice = slice(self.ipntr[0]-1, self.ipntr[0]-1+self.n)
+        yslice = slice(self.ipntr[1]-1, self.ipntr[1]-1+self.n)
+        if self.ido == -1:
+            # initialization
+            self.workd[yslice] = self.matvec(self.workd[xslice])
+        elif self.ido == 1:
+            # compute y=Ax
+            self.workd[yslice] = self.matvec(self.workd[xslice])
+        else:
+            self.converged = True
+
+            if self.info < -1 :
+                raise RuntimeError("Error info=%d in arpack" % self.info)
+            elif self.info == -1:
+                warnings.warn("Maximum number of iterations taken: %s" % self.iparam[2])
+
+            if self.iparam[4] < self.k:
+                warnings.warn("Only %d/%d eigenvectors converged" % (self.iparam[4], self.k))
+
 def eigen(A, k=6, M=None, sigma=None, which='LM', v0=None,
           ncv=None, maxiter=None, tol=0,
           return_eigenvectors=True):
@@ -217,42 +276,16 @@
         raise ValueError('expected square matrix (shape=%s)' % A.shape)
     n = A.shape[0]
 
-    params = _ArpackParams(n, k, A.dtype.char, "unsymmetric", sigma,
+    matvec = lambda x : A.matvec(x)
+    params = _ArpackParams(n, k, A.dtype.char, matvec, "unsymmetric", sigma,
                            ncv, v0, maxiter, which, tol)
 
     if M is not None:
         raise NotImplementedError("generalized eigenproblem not supported yet")
 
-    ido = 0
+    while not params.converged:
+        params.iterate()
 
-    while True:
-        if params.tp in 'fd':
-            ido, params.resid, params.v, params.iparam, params.ipntr, params.info = \
-                params.solver(ido, params.bmat, params.which, params.k, params.tol,
-                        params.resid, params.v, params.iparam, params.ipntr,
-                        params.workd, params.workl, params.info)
-        else:
-            ido, params.resid, params.v, params.iparam, params.ipntr, params.info =\
-                params.solver(ido, params.bmat, params.which, params.k, params.tol,
-                        params.resid, params.v, params.iparam, params.ipntr,
-                        params.workd, params.workl, params.rwork, params.info)
-
-        xslice = slice(params.ipntr[0]-1, params.ipntr[0]-1+n)
-        yslice = slice(params.ipntr[1]-1, params.ipntr[1]-1+n)
-        if ido == -1:
-            # initialization
-            params.workd[yslice] = A.matvec(params.workd[xslice])
-        elif ido == 1:
-            # compute y=Ax
-            params.workd[yslice] = A.matvec(params.workd[xslice])
-        else:
-            break
-
-    if params.info < -1 :
-        raise RuntimeError("Error info=%d in arpack" % params.info)
-    elif params.info == -1:
-        warnings.warn("Maximum number of iterations taken: %s" % self.iparam[2])
-
     # now extract eigenvalues and (optionally) eigenvectors
     rvec = return_eigenvectors
     ierr = 0
@@ -419,35 +452,13 @@
     if M is not None:
         raise NotImplementedError("generalized eigenproblem not supported yet")
 
-    params = _ArpackParams(n, k, A.dtype.char, "symmetric", sigma,
+    matvec = lambda x : A.matvec(x)
+    params = _ArpackParams(n, k, A.dtype.char, matvec, "symmetric", sigma,
                            ncv, v0, maxiter, which, tol)
 
-    ido = 0
-    while True:
-        ido, params.resid, params.v, params.iparam, params.ipntr, params.info = \
-            params.solver(ido, params.bmat, params.which, params.k, params.tol,
-                    params.resid, params.v, params.iparam, params.ipntr,
-                    params.workd, params.workl, params.info)
+    while not params.converged:
+        params.iterate()
 
-        xslice = slice(params.ipntr[0]-1, params.ipntr[0]-1+n)
-        yslice = slice(params.ipntr[1]-1, params.ipntr[1]-1+n)
-        if ido == -1:
-            # initialization
-            params.workd[yslice] = A.matvec(params.workd[xslice])
-        elif ido == 1:
-            # compute y=Ax
-            params.workd[yslice] = A.matvec(params.workd[xslice])
-        else:
-            break
-
-    if params.info < -1 :
-        raise RuntimeError("Error info=%d in arpack" % params.info)
-    elif params.info == 1:
-        warnings.warn("Maximum number of iterations taken: %s" % params.iparam[2])
-
-    if params.iparam[4] < k:
-        warnings.warn("Only %d/%d eigenvectors converged" % (params.iparam[4], k))
-
     # now extract eigenvalues and (optionally) eigenvectors
     rvec = return_eigenvectors
     ierr = 0



More information about the Scipy-svn mailing list