[Scipy-svn] r6987 - in trunk/scipy/sparse/linalg/eigen/arpack: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Dec 4 14:25:04 CST 2010


Author: ptvirtan
Date: 2010-12-04 14:25:03 -0600 (Sat, 04 Dec 2010)
New Revision: 6987

Modified:
   trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
   trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py
Log:
BUG: sparse/arpack: fix ARPACK error handling

Make exceptions instances of ArpackError. Raise ArpackNoConvergence
always when no convergence is obtained.

Modified: trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-12-04 20:24:49 UTC (rev 6986)
+++ trunk/scipy/sparse/linalg/eigen/arpack/arpack.py	2010-12-04 20:25:03 UTC (rev 6987)
@@ -39,10 +39,8 @@
 
 __docformat__ = "restructuredtext en"
 
-__all___=['eigen','eigen_symmetric', 'svd']
+__all___=['eigen','eigen_symmetric', 'svd', 'ArpackNoConvergence']
 
-import warnings
-
 import _arpack
 import numpy as np
 from scipy.sparse.linalg.interface import aslinearoperator
@@ -51,13 +49,123 @@
 _type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
 _ndigits = {'f':5, 'd':12, 'F':5, 'D':12}
 
+_NAUPD_ERRORS = {
+    0: "Normal exit.",
+    1: "Maximum number of iterations taken. "
+       "All possible eigenvalues of OP has been found.",
+    2: "No longer an informational error. Deprecated starting with "
+       "release 2 of ARPACK.",
+    3: "No shifts could be applied during a cycle of the Implicitly "
+       "restarted Arnoldi iteration. One possibility is to increase "
+       "the size of NCV relative to NEV. ",
+    -1: "N must be positive.",
+    -2: "NEV must be positive.",
+    -3: "NCV must be greater than NEV and less than or equal to N.",
+    -4: "The maximum number of Arnoldi update iterations allowed "
+        "must be greater than zero.",
+    -5: "WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.",
+    -6: "BMAT must be one of 'I' or 'G'.",
+    -7: "Length of private work array WORKL is not sufficient.",
+    -8: "Error return from trid. eigenvalue calculation; "
+        "Informational error from LAPACK routine dsteqr .",
+    -9: "Starting vector is zero.",
+    -10: "IPARAM(7) must be 1,2,3,4,5.",
+    -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatable.",
+    -12: "IPARAM(1) must be equal to 0 or 1.",
+    -13: "NEV and WHICH = 'BE' are incompatable. ",
+    -9999: "Could not build an Arnoldi factorization. "
+           "IPARAM(5) returns the size of the current Arnoldi "
+           "factorization. The user is advised to check that "
+           "enough workspace and array storage has been allocated.",
+}
+
+_NEUPD_ERRORS = {
+    0: "Normal exit.",
+    1: "The Schur form computed by LAPACK routine dlahqr "
+       "could not be reordered by LAPACK routine dtrsen. "
+       "Re-enter subroutine dneupd  with IPARAM(5)NCV and "
+       "increase the size of the arrays DR and DI to have "
+       "dimension at least dimension NCV and allocate at least NCV "
+       "columns for Z. NOTE: Not necessary if Z and V share "
+       "the same space. Please notify the authors if this error"
+       "occurs.",
+    -1: "N must be positive.",
+    -2: "NEV must be positive.",
+    -3: "NCV-NEV >= 2 and less than or equal to N.",
+    -5: "WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'",
+    -6: "BMAT must be one of 'I' or 'G'.",
+    -7: "Length of private work WORKL array is not sufficient.",
+    -8: "Error return from calculation of a real Schur form. "
+        "Informational error from LAPACK routine dlahqr .",
+    -9: "Error return from calculation of eigenvectors. "
+        "Informational error from LAPACK routine dtrevc.",
+    -10: "IPARAM(7) must be 1,2,3,4.",
+    -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.",
+    -12: "HOWMNY = 'S' not yet implemented",
+    -13: "HOWMNY must be one of 'A' or 'P' if RVEC = .true.",
+    -14: "DNAUPD  did not find any eigenvalues to sufficient "
+         "accuracy.",
+    -15: "DNEUPD got a different count of the number of converged "
+         "Ritz values than DNAUPD got.  This indicates the user "
+         "probably made an error in passing data from DNAUPD to "
+         "DNEUPD or that the data was modified before entering "
+         "DNEUPD",
+}
+
+_SEUPD_ERRORS = {
+    0: "Normal exit.",
+    -1: "N must be positive.",
+    -2: "NEV must be positive.",
+    -3: "NCV must be greater than NEV and less than or equal to N.",
+    -5: "WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.",
+    -6: "BMAT must be one of 'I' or 'G'.",
+    -7: "Length of private work WORKL array is not sufficient.",
+    -8: ("Error return from trid. eigenvalue calculation; "
+         "Information error from LAPACK routine dsteqr."),
+    -9: "Starting vector is zero.",
+    -10: "IPARAM(7) must be 1,2,3,4,5.",
+    -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.",
+    -12: "NEV and WHICH = 'BE' are incompatible.",
+    -14: "DSAUPD  did not find any eigenvalues to sufficient accuracy.",
+    -15: "HOWMNY must be one of 'A' or 'S' if RVEC = .true.",
+    -16: "HOWMNY = 'S' not yet implemented",
+    -17: ("DSEUPD  got a different count of the number of converged "
+          "Ritz values than DSAUPD  got.  This indicates the user "
+          "probably made an error in passing data from DSAUPD  to "
+          "DSEUPD  or that the data was modified before entering  "
+          "DSEUPD.")
+}
+
+class ArpackError(RuntimeError):
+    """
+    ARPACK error
+    """
+    def __init__(self, info, infodict=_NAUPD_ERRORS):
+        msg = infodict.get(info, "Unknown error")
+        RuntimeError.__init__(self, "ARPACK error %d: %s" % (info, msg))
+
+class ArpackNoConvergence(ArpackError):
+    """
+    ARPACK iteration did not converge
+
+    Attributes
+    ----------
+    eigenvalues : ndarray
+        Partial result. Converged eigenvalues.
+    eigenvectors : ndarray
+        Partial result. Converged eigenvectors.
+
+    """
+    def __init__(self, msg, eigenvalues, eigenvectors):
+        ArpackError.__init__(self, -1, {-1: msg})
+        self.eigenvalues = eigenvalues
+        self.eigenvectors = eigenvectors
+
 class _ArpackParams(object):
     def __init__(self, n, k, tp, matvec, sigma=None,
                  ncv=None, v0=None, maxiter=None, which="LM", tol=0):
         if k <= 0:
             raise ValueError("k must be positive, k=%d" % k)
-        if k == n:
-            raise ValueError("k must be less than rank(A), k=%d" % k)
 
         if maxiter is None:
             maxiter = n * 10
@@ -107,11 +215,26 @@
         self.converged = False
         self.ido = 0
 
+    def _raise_no_convergence(self):
+        msg = "No convergence (%d iterations, %d/%d eigenvectors converged)"
+        k_ok = self.iparam[4]
+        num_iter = self.iparam[2]
+        try:
+            ev, vec = self.extract(True)
+        except ArpackError, err:
+            msg = "%s [%s]" % (msg, err)
+            ev = np.zeros((0,))
+            vec = np.zeros((self.n, 0))
+            k_ok = 0
+        raise ArpackNoConvergence(msg % (num_iter, k_ok, self.k), ev, vec)
+
 class _SymmetricArpackParams(_ArpackParams):
     def __init__(self, n, k, tp, matvec, sigma=None,
                  ncv=None, v0=None, maxiter=None, which="LM", tol=0):
         if not which in ['LM', 'SM', 'LA', 'SA', 'BE']:
             raise ValueError("which must be one of %s" % ' '.join(whiches))
+        if k >= n:
+            raise ValueError("k must be less than rank(A), k=%d" % k)
 
         _ArpackParams.__init__(self, n, k, tp, matvec, sigma,
                  ncv, v0, maxiter, which, tol)
@@ -145,14 +268,13 @@
         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.info == 0:
+                pass
+            elif self.info == 1:
+                self._raise_no_convergence()
+            else:
+                raise ArpackError(self.info)
 
-            if self.iparam[4] < self.k:
-                warnings.warn("Only %d/%d eigenvectors converged" % (self.iparam[4], self.k))
-
     def extract(self, return_eigenvectors):
         rvec = return_eigenvectors
         ierr = 0
@@ -160,14 +282,18 @@
         sselect = np.zeros(self.ncv, 'int') # unused
         sigma = 0.0 # no shifts, not implemented
 
-        d, z, info = self._arpack_extract(rvec, howmny, sselect, sigma, self.bmat,
+        d, z, ierr = self._arpack_extract(rvec, howmny, sselect, sigma, self.bmat,
                 self.which, self.k, self.tol, self.resid, self.v,
                 self.iparam[0:7], self.ipntr, self.workd[0:2*self.n],
                 self.workl,ierr)
 
         if ierr != 0:
-            raise RuntimeError("Error info=%d in arpack" % params.info)
+            raise ArpackError(ierr, infodict=_SEUPD_ERRORS)
 
+        k_ok = self.iparam[4]
+        d = d[:k_ok]
+        z = z[:,:k_ok]
+
         if return_eigenvectors:
             return d, z
         else:
@@ -178,6 +304,8 @@
                  ncv=None, v0=None, maxiter=None, which="LM", tol=0):
         if not which in ["LM", "SM", "LR", "SR", "LI", "SI"]:
             raise ValueError("Parameter which must be one of %s" % ' '.join(whiches))
+        if k >= n-1:
+            raise ValueError("k must be less than rank(A)-1, k=%d" % k)
 
         _ArpackParams.__init__(self, n, k, tp, matvec, sigma,
                  ncv, v0, maxiter, which, tol)
@@ -222,10 +350,12 @@
         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.info == 0:
+                pass
+            elif self.info == 1:
+                self._raise_no_convergence()
+            else:
+                raise ArpackError(self.info)
 
     def extract(self, return_eigenvectors):
         k, n = self.k, self.n
@@ -241,13 +371,16 @@
             dr = np.zeros(k+1, self.tp)
             di = np.zeros(k+1, self.tp)
             zr = np.zeros((n, k+1), self.tp)
-            dr, di, zr, self.info=\
+            dr, di, zr, ierr=\
                 self._arpack_extract(return_eigenvectors,
                        howmny, sselect, sigmar, sigmai, workev,
                        self.bmat, self.which, k, self.tol, self.resid,
                        self.v, self.iparam, self.ipntr,
                        self.workd, self.workl, self.info)
 
+            if ierr != 0:
+                raise ArpackError(ierr, infodict=_NEUPD_ERRORS)
+
             # The ARPACK nonsymmetric real and double interface (s,d)naupd return
             # eigenvalues and eigenvectors in real (float,double) arrays.
 
@@ -294,16 +427,21 @@
 
         else:
             # complex is so much simpler...
-            d, z, self.info =\
+            d, z, ierr =\
                     self._arpack_extract(return_eigenvectors,
                            howmny, sselect, sigmar, workev,
                            self.bmat, self.which, k, self.tol, self.resid,
                            self.v, self.iparam, self.ipntr,
                            self.workd, self.workl, self.rwork, ierr)
 
-        if ierr != 0:
-            raise RuntimeError("Error info=%d in arpack" % info)
+            if ierr != 0:
+                raise ArpackError(ierr, infodict=_NEUPD_ERRORS)
 
+            k_ok = self.iparam[4]
+            d = d[:k_ok]
+            z = z[:,:k_ok]
+
+
         if return_eigenvectors:
             return d, z
         else:
@@ -325,7 +463,9 @@
         the matrix vector product A * x.  The sparse matrix formats
         in scipy.sparse are appropriate for A.
     k : integer
-        The number of eigenvalues and eigenvectors desired
+        The number of eigenvalues and eigenvectors desired.
+        `k` must be smaller than N. It is not possible to compute all
+        eigenvectors of a matrix.
 
     Returns
     -------
@@ -363,9 +503,19 @@
         Maximum number of Arnoldi update iterations allowed
     tol : float
         Relative accuracy for eigenvalues (stopping criterion)
+        The default value of 0 implies machine precision.
     return_eigenvectors : boolean
         Return eigenvectors (True) in addition to eigenvalues
 
+    Raises
+    ------
+    ArpackNoConvergence
+        When the requested convergence is obtained.
+
+        The currently converged eigenvalues and eigenvectors can be found
+        as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
+        objed.
+
     See Also
     --------
     eigen_symmetric : eigenvalues and eigenvectors for symmetric matrix A
@@ -445,7 +595,8 @@
 
     ncv : integer
         The number of Lanczos vectors generated
-        ncv must be greater than k; it is recommended that ncv > 2*k
+        ncv must be greater than k and smaller than n;
+        it is recommended that ncv > 2*k
 
     which : string
         Which k eigenvectors and eigenvalues to find:
@@ -460,11 +611,21 @@
         Maximum number of Arnoldi update iterations allowed
 
     tol : float
-        Relative accuracy for eigenvalues (stopping criterion)
+        Relative accuracy for eigenvalues (stopping criterion).
+        The default value of 0 implies machine precision.
 
     return_eigenvectors : boolean
         Return eigenvectors (True) in addition to eigenvalues
 
+    Raises
+    ------
+    ArpackNoConvergence
+        When the requested convergence is obtained.
+
+        The currently converged eigenvalues and eigenvectors can be found
+        as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
+        objed.
+
     See Also
     --------
     eigen : eigenvalues and eigenvectors for a general (nonsymmetric) matrix A

Modified: trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py	2010-12-04 20:24:49 UTC (rev 6986)
+++ trunk/scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py	2010-12-04 20:25:03 UTC (rev 6987)
@@ -11,8 +11,9 @@
         assert_raises, verbose
 
 from numpy import array, finfo, argsort, dot, round, conj, random
-from scipy.sparse import csc_matrix
-from scipy.sparse.linalg.eigen.arpack import eigen_symmetric, eigen, svd
+from scipy.sparse import csc_matrix, isspmatrix
+from scipy.sparse.linalg.eigen.arpack import eigen_symmetric, eigen, svd, \
+     ArpackNoConvergence
 
 from scipy.linalg import svd as dsvd
 
@@ -132,6 +133,23 @@
             self.eval_evec(self.symmetric[0],typ,k,which='LM',v0=v0)
 
 
+    def test_no_convergence(self):
+        np.random.seed(1234)
+        m = np.random.rand(30, 30)
+        m = m + m.T
+        try:
+            w, v = eigen_symmetric(m, 4, which='LM', v0=m[:,0], maxiter=5)
+            raise AssertionError("Spurious no-error exit")
+        except ArpackNoConvergence, err:
+            k = len(err.eigenvalues)
+            if k <= 0:
+                raise AssertionError("Spurious no-eigenvalues-found case")
+            w, v = err.eigenvalues, err.eigenvectors
+            for ww, vv in zip(w, v.T):
+                assert_array_almost_equal(dot(m, vv), ww*vv,
+                                          decimal=_ndigits['d'])
+
+
 class TestEigenComplexSymmetric(TestArpack):
 
     def sort_choose(self,eval,typ,k,which):
@@ -173,6 +191,21 @@
                 self.eval_evec(self.symmetric[0],typ,k,which)
 
 
+    def test_no_convergence(self):
+        np.random.seed(1234)
+        m = np.random.rand(30, 30) + 1j*np.random.rand(30, 30)
+        try:
+            w, v = eigen(m, 3, which='LM', v0=m[:,0], maxiter=30)
+            raise AssertionError("Spurious no-error exit")
+        except ArpackNoConvergence, err:
+            k = len(err.eigenvalues)
+            if k <= 0:
+                raise AssertionError("Spurious no-eigenvalues-found case")
+            w, v = err.eigenvalues, err.eigenvectors
+            for ww, vv in zip(w, v.T):
+                assert_array_almost_equal(dot(m, vv), ww*vv,
+                                          decimal=_ndigits['D'])
+
 class TestEigenNonSymmetric(TestArpack):
 
 
@@ -231,9 +264,21 @@
             v0 = random.rand(n).astype(typ)
             self.eval_evec(self.symmetric[0],typ,k,which='LM',v0=v0)
 
+    def test_no_convergence(self):
+        np.random.seed(1234)
+        m = np.random.rand(30, 30)
+        try:
+            w, v = eigen(m, 3, which='LM', v0=m[:,0], maxiter=30)
+            raise AssertionError("Spurious no-error exit")
+        except ArpackNoConvergence, err:
+            k = len(err.eigenvalues)
+            if k <= 0:
+                raise AssertionError("Spurious no-eigenvalues-found case")
+            w, v = err.eigenvalues, err.eigenvectors
+            for ww, vv in zip(w, v.T):
+                assert_array_almost_equal(dot(m, vv), ww*vv,
+                                          decimal=_ndigits['d'])
 
-
-
 class TestEigenComplexNonSymmetric(TestArpack):
 
     def sort_choose(self,eval,typ,k,which):
@@ -287,6 +332,21 @@
                     self.eval_evec(m,typ,k,which)
 
 
+    def test_no_convergence(self):
+        np.random.seed(1234)
+        m = np.random.rand(30, 30) + 1j*np.random.rand(30, 30)
+        try:
+            w, v = eigen(m, 3, which='LM', v0=m[:,0], maxiter=30)
+            raise AssertionError("Spurious no-error exit")
+        except ArpackNoConvergence, err:
+            k = len(err.eigenvalues)
+            if k <= 0:
+                raise AssertionError("Spurious no-eigenvalues-found case")
+            w, v = err.eigenvalues, err.eigenvectors
+            for ww, vv in zip(w, v.T):
+                assert_array_almost_equal(dot(m, vv), ww*vv,
+                                          decimal=_ndigits['D'])
+
 def test_eigen_bad_shapes():
     # A is not square.
     A = csc_matrix(np.zeros((2,3)))



More information about the Scipy-svn mailing list