[Scipy-svn] r6987 - in trunk/scipy/sparse/linalg/eigen/arpack: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Dec 4 15:25:04 EST 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