# [Scipy-svn] r2324 - in trunk/Lib/sandbox/arpack: . tests

scipy-svn at scipy.org scipy-svn at scipy.org
Mon Nov 20 14:29:04 CST 2006

```Author: nmarais
Date: 2006-11-20 14:27:36 -0600 (Mon, 20 Nov 2006)
New Revision: 2324

Modified:
trunk/Lib/sandbox/arpack/speigs.py
trunk/Lib/sandbox/arpack/tests/test_speigs.py
Log:
Refactor lower-level ARPACK drivers and improve docstrings

Modified: trunk/Lib/sandbox/arpack/speigs.py
===================================================================
--- trunk/Lib/sandbox/arpack/speigs.py	2006-11-18 02:09:53 UTC (rev 2323)
+++ trunk/Lib/sandbox/arpack/speigs.py	2006-11-20 20:27:36 UTC (rev 2324)
@@ -1,49 +1,56 @@
import numpy as N
import _arpack
+import warnings

-def eigvals(matvec, n, nev, ncv=None):
-    """
-    Calculate eigenvalues for system with matrix-vector product matvec, dimension n
+__all___=['ArpackException','ARPACK_eigs', 'ARPACK_gen_eigs']

-    Arguments
-    =========
-    matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> A*x
-    n -- Matrix dimension of the problem
-    nev -- Number of eigenvalues to calculate
-    ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev
-
-    Return Values
-    =============
-    Real array of nev eigenvalues, or complex array if values are complex
-    """
-
+class ArpackException(RuntimeError):
+    ARPACKErrors = { 0: """Normal exit.""",
+                     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-NEV >= 2 and less than or equal to N.""",
+                     -4: """The maximum number of Arnoldi update iteration
+                     must be greater than zero.""",
+                     -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 array is not sufficient.""",
+                     -8: """Error return from LAPACK eigenvalue calculation;""",
+                     -9: """Starting vector is zero.""",
+                     -10: """IPARAM(7) must be 1,2,3,4.""",
+                     -11: """IPARAM(7) = 1 and BMAT = 'G' are incompatable.""",
+                     -12: """IPARAM(1) must be equal to 0 or 1.""",
+                     -9999: """Could not build an Arnoldi factorization.
+                     IPARAM(5) returns the size of the current Arnoldi
+                     factorization.""",
+                     }
+    def __init__(self, info):
+        self.info = info
+    def __str__(self):
+        try: return self.ARPACKErrors[self.info]
+        except KeyError: return "Unknown ARPACK error"
+
+def check_init(n, nev, ncv):
assert(nev <= n-4)  # ARPACK seems to cause a segfault otherwise
if ncv is None:
-        ncv = min(2*nev, n-1)
-    iparam = N.zeros(11, N.int32) # Array with assorted extra paramters for F77 call
-    ishfts = 1         # Some random arpack parameter
-    maxitr = n*3       # Maximum number of iterations
-    # Some random arpack parameter (I think it tells ARPACK to solve the
-    # regular eigenproblem the "standard" way
-    mode1 = 1
-    iparam[0] = ishfts
-    iparam[2] = maxitr
-    iparam[6] = mode1
+        ncv = min(2*nev+1, n-1)
+    maxitr = max(n, 1000)       # Maximum number of iterations
+    return ncv, maxitr
+
+def init_workspaces(n,nev,ncv):
ipntr = N.zeros(14, N.int32) # Pointers into memory structure used by F77 calls
d = N.zeros((ncv, 3), N.float64, order='FORTRAN') # Temp workspace
# Temp workspace/error residuals upon iteration completion
resid = N.zeros(n, N.float64)
-
workd = N.zeros(3*n, N.float64) # workspace
workl = N.zeros(3*ncv*ncv+6*ncv, N.float64) # more workspace
-    tol = 1e-16            # Desired convergence tollerance
-    ido = 0                # Communication variable used by ARPACK to tell the user what to do
-    info = 0               # Used for error reporting
-    which = 'SM'           # Request smallest magnitude eigenvalues, see dnaupd.f for more options
-    bmat  = 'I'            # Standard (not generalised) eigenproblem
# Storage for the Arnoldi basis vectors
-    v = N.zeros((n, ncv), dtype=float, order='FORTRAN')
-
+    v = N.zeros((n, ncv), dtype=N.float64, order='FORTRAN')
+    return (ipntr, d, resid, workd, workl, v)
+
+def init_debug():
# Causes various debug info to be printed by ARPACK
_arpack.debug.ndigit = -3
_arpack.debug.logfil = 6
@@ -54,43 +61,70 @@
_arpack.debug.mneigh = 0
_arpack.debug.mneupd = 1

-
-    cnt = 0
-    # Arnouldi iteration.
-    while True:
-        ido,resid,v,iparam,ipntr,info = _arpack.dnaupd(
-            ido, bmat, which, nev, tol, resid, v, iparam, ipntr, workd, workl, info)
-        # Exit if reverse communication flag does not request a matrix-vector
-        # product
-        if ido not in  (-1, 1): break
-        cnt += 1
-        x = workd[ipntr[0]-1:ipntr[0]+n-1]
-        workd[ipntr[1]-1:ipntr[1]+n-1] = matvec(x)    # y = A*x
-
-
-    if info != 0: raise "Hell" # Indicates some error during the Arnouldi iterations
-
-    dr = N.zeros(nev+1, N.float64) # Real part of eigenvalues
-    di = N.zeros(nev+1, N.float64) # Imaginary part of eigenvalues
+def init_postproc_workspace(n, nev, ncv):
# Used as workspace and to return eigenvectors if requested. Not touched if
# eigenvectors are not requested
-    z = N.zeros((n, nev+1), N.float64, order='FORTRAN')
workev = N.zeros(3*ncv, N.float64) # More workspace
select = N.zeros(ncv, N.int32) # Used as internal workspace since dneupd
# parameter HOWMNY == 'A'
+    return (workev, select)

+def postproc(n, nev, ncv, sigmar, sigmai, bmat, which,
+             tol, resid, v, iparam, ipntr, workd, workl, info):
+    workev, select = init_postproc_workspace(n, nev, ncv)
+    ierr = 0
# Postprocess the Arnouldi vectors to extract eigenvalues/vectors
# If dneupd's first paramter is 'True' the eigenvectors are also calculated,
# 'False' only the eigenvalues
dr,di,z,info = _arpack.dneupd(
-        True, 'A', select, 0., 0., workev, bmat, which, nev, tol,
-        resid, v, iparam, ipntr, workd, workl, info)
+        True, 'A', select, sigmar, sigmai, workev, bmat, which, nev, tol, resid, v,
+        iparam, ipntr, workd, workl, info)
+
if N.abs(di[:-1]).max() == 0: dr = dr[:-1]
else: dr =  dr[:-1] + 1j*di[:-1]
return (dr, z[:,:-1])

-def geneigvals(matvec, sigma_solve, n, sigma, nev, ncv=None):
+
+def ARPACK_eigs(matvec, n, nev, which='SM', ncv=None, tol=1e-14):
"""
+    Calculate eigenvalues for system with matrix-vector product matvec, dimension n
+
+    Arguments
+    =========
+    matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> A*x
+    n -- Matrix dimension of the problem
+    nev -- Number of eigenvalues to calculate
+    which -- Spectrum selection. See details below. Defaults to 'SM'
+    ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev+1
+    tol -- Numerical tollerance for Arnouldi iteration convergence. Defaults to 1e-14
+
+    Spectrum Selection
+    ==================
+    which can take one of several values:
+
+    'LM' -> Request eigenvalues with largest magnitude.
+    'SM' -> Request eigenvalues with smallest magnitude.
+    'LR' -> Request eigenvalues with largest real part.
+    'SR' -> Request eigenvalues with smallest real part.
+    'LI' -> Request eigenvalues with largest imaginary part.
+    'SI' -> Request eigenvalues with smallest imaginary part.
+
+    Return Values
+    =============
+    (eig_vals, eig_vecs) where eig_vals are the requested eigenvalues and
+    eig_vecs the corresponding eigenvectors. If all the eigenvalues are real,
+    eig_vals is a real array but if some eigenvalues are complex it is a
+    complex array.
+
+    """
+    bmat = 'I'                          # Standard eigenproblem
+    ncv, resid, iparam, ipntr, v, workd, workl, info = ARPACK_iteration(
+        matvec, lambda x: x, n, bmat, which, nev, tol, ncv, mode=1)
+    return postproc(n, nev, ncv, 0., 0., bmat, which, tol,
+                    resid, v, iparam, ipntr, workd, workl, info)
+
+def ARPACK_gen_eigs(matvec, sigma_solve, n, sigma, nev, which='LR', ncv=None, tol=1e-14):
+    """
Calculate eigenvalues close to sigma for generalised eigen system

Given a system [A]x = k_i*[M]x where [A] and [M] are matrices and k_i are
@@ -104,53 +138,73 @@
n -- Matrix dimension of the problem
sigma -- Eigenvalue spectral shift real value
nev -- Number of eigenvalues to calculate
-    ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev
+    which -- Spectrum selection. See details below. Defaults to 'LR'
+    ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev+1
+    tol -- Numerical tollerance for Arnouldi iteration convergence. Defaults to 1e-14

+    Spectrum Shift
+    ==============
+
+    The spectrum of the orignal system is shifted by sigma. This transforms the
+    original eigenvalues to be 1/(original_eig-sigma) in the shifted
+    system. ARPACK then operates on the shifted system, transforming it back to
+    the original system in a postprocessing step.
+
+    The spectrum shift causes eigenvalues close to sigma to become very large
+    in the transformed system. This allows quick convergence for these
+    eigenvalues. This is particularly useful if a system has a number of
+    trivial zero-eigenvalues that are to be ignored.
+
+    Spectrum Selection
+    ==================
+    which can take one of several values:
+
+    'LM' -> Request spectrum shifted eigenvalues with largest magnitude.
+    'SM' -> Request spectrum shifted eigenvalues with smallest magnitude.
+    'LR' -> Request spectrum shifted eigenvalues with largest real part.
+    'SR' -> Request spectrum shifted eigenvalues with smallest real part.
+    'LI' -> Request spectrum shifted eigenvalues with largest imaginary part.
+    'SI' -> Request spectrum shifted eigenvalues with smallest imaginary part.
+
+    The effect on the actual system is:
+    'LM' -> Eigenvalues closest to sigma on the complex plane
+    'LR' -> Eigenvalues with real part > sigma, provided they exist
+
+
Return Values
=============
-    Real array of nev eigenvalues, or complex array if values are complex
+    (eig_vals, eig_vecs) where eig_vals are the requested eigenvalues and
+    eig_vecs the corresponding eigenvectors. If all the eigenvalues are real,
+    eig_vals is a real array but if some eigenvalues are complex it is a
+    complex array. The eigenvalues and vectors correspond to the original
+    system, not the shifted system. The shifted system is only used interally.
+
"""
-    if ncv is None:
-        ncv = min(2*nev, n-1)
+    bmat = 'G'                          # Generalised eigenproblem
+    ncv, resid, iparam, ipntr, v, workd, workl, info = ARPACK_iteration(
+        matvec, sigma_solve, n, bmat, which, nev, tol, ncv, mode=3)
+    sigmar = sigma
+    sigmai = 0.
+    return postproc(n, nev, ncv, sigmar, sigmai, bmat, which, tol,
+                    resid, v, iparam, ipntr, workd, workl, info)

-    iparam = N.zeros(11, N.int32) # Array with assorted extra paramters for F77 call
+def ARPACK_iteration(matvec, sigma_solve, n, bmat, which, nev, tol, ncv, mode):
+    ncv, maxitr = check_init(n, nev, ncv)
+    ipntr, d, resid, workd, workl, v = init_workspaces(n,nev,ncv)
+    init_debug()
ishfts = 1         # Some random arpack parameter
-    maxitr = n*3       # Maximum number of iterations
# Some random arpack parameter (I think it tells ARPACK to solve the
-    # regular eigenproblem the "standard" way
-    mode = 3
-    iparam[0] = ishfts
-    iparam[2] = maxitr
-    iparam[6] = mode
-    ipntr = N.zeros(14, N.int32) # Pointers into memory structure used by F77 calls
-    d = N.zeros((ncv, 3), N.float64, order='FORTRAN') # Temp workspace
-    # Temp workspace/error residuals upon iteration completion
-    resid = N.zeros(n, N.float64)
-
-    workd = N.zeros(3*n, N.float64) # workspace
-    workl = N.zeros(3*ncv*ncv+6*ncv, N.float64) # more workspace
-    tol = 1e-7             # Desired convergence tollerance
+    # general eigenproblem using shift-invert
+    iparam = N.zeros(11, N.int32) # Array with assorted extra paramters for F77 call
+    iparam[[0,2,6]] = ishfts, maxitr, mode
ido = 0                # Communication variable used by ARPACK to tell the user what to do
info = 0               # Used for error reporting
-    which = 'LR'           # Request largest magnitude eigenvalues, see dnaupd.f for more options
-    bmat  = 'G'            # Generalised eigenproblem
-    # Storage for the Arnoldi basis vectors
-    v = N.zeros((n, ncv), dtype=float, order='FORTRAN')
-    # Causes various debug info to be printed by ARPACK
-    _arpack.debug.ndigit = -3
-    _arpack.debug.logfil = 6
-    _arpack.debug.mnaitr = 0
-    _arpack.debug.mnapps = 0
-    _arpack.debug.mnaupd = 1
-    _arpack.debug.mnaup2 = 0
-    _arpack.debug.mneigh = 0
-    _arpack.debug.mneupd = 1
-
# Arnouldi iteration.
while True:
-        ido,resid,v,iparam,ipntr,info = _arpack.dnaupd(ido, bmat, which, nev, tol, resid, v,
-                            iparam, ipntr, workd, workl, info)
-        if ido == -1:  # Perform y = inv[A - sigma*M]*M*x
+        ido,resid,v,iparam,ipntr,info = _arpack.dnaupd(
+            ido, bmat, which, nev, tol, resid, v, iparam, ipntr, workd, workl, info)
+        if ido == -1 or ido == 1 and mode not in (3,4):
+            # Perform y = inv[A - sigma*M]*M*x
x = workd[ipntr[0]-1:ipntr[0]+n-1]
Mx = matvec(x)    # Mx = [M]*x
workd[ipntr[1]-1:ipntr[1]+n-1] = sigma_solve(Mx)
@@ -163,26 +217,9 @@
workd[ipntr[1]-1:ipntr[1]+n-1] = matvec(x)
else:          # Finished, or error
break
-        if info != 0: raise "Hell" # Indicates some error during the Arnouldi iterations
+        if info == 1:
+            warn.warn("Maximum number of iterations taken: %s"%iparam[2])
+        elif info != 0:
+            raise ArpackException(info)

-    # Used as workspace and to return eigenvectors if requested. Not touched if
-    # eigenvectors are not requested
-    z = N.zeros((n, nev+1), N.float64, order='FORTRAN')
-    workev = N.zeros(3*ncv, N.float64) # More workspace
-    ierr = 0
-    select = N.zeros(ncv, N.int32) # Used as internal workspace since dneupd
-                                   # parameter HOWMNY == 'A'
-    sigmar = sigma
-    sigmai = 0.
-    # Postprocess the Arnouldi vectors to extract eigenvalues/vectors
-    # If dneupd's first paramter is 'True' the eigenvectors are also calculated,
-    # 'False' only the eigenvalues
-    dr,di,z,info = _arpack.dneupd(
-        True, 'A', select, sigmar, sigmai, workev, bmat, which, nev, tol, resid, v,
-        iparam, ipntr, workd, workl, info)
-
-    if N.abs(di[:-1]).max() == 0:
-        return dr[:-1]
-    else:
-        return dr[:-1] + 1j*di[:-1]
-
+    return (ncv, resid, iparam, ipntr, v, workd, workl, info)

Modified: trunk/Lib/sandbox/arpack/tests/test_speigs.py
===================================================================
--- trunk/Lib/sandbox/arpack/tests/test_speigs.py	2006-11-18 02:09:53 UTC (rev 2323)
+++ trunk/Lib/sandbox/arpack/tests/test_speigs.py	2006-11-20 20:27:36 UTC (rev 2324)
@@ -28,7 +28,7 @@
matvec = get_matvec(A)
#= lambda x: N.asarray(A*x)[0]
nev=4
-        eigvs = eigvals(matvec, A.shape[0], nev=nev)
+        eigvs = ARPACK_eigs(matvec, A.shape[0], nev=nev)
calc_vals = eigvs[0]
# Ensure the calculate eigenvectors have the same sign as the refence values
calc_vecs = eigvs[1] / [N.sign(x[0]) for x in eigvs[1].T]
@@ -40,11 +40,14 @@
#     def test(self):
#         import pickle
#         import scipy.linsolve
#         sigma = 27.
#         sigma_solve = scipy.linsolve.splu(A - sigma*B).solve
-#         w = geneigvals(B.matvec, sigma_solve, B.shape[0], sigma, 10)
-
+#         w = ARPACK_gen_eigs(B.matvec, sigma_solve, B.shape[0], sigma, 10)[0]
+#         assert_array_almost_equal(w,
+#         [27.346442255386375,  49.100299170945405,  56.508474856551544, 56.835800191692492,
+#          65.944215785041365, 66.194792400328367, 78.003788872725238, 79.550811647295944,
+#          94.646308846854879, 95.30841709116271], decimal=11)

if __name__ == "__main__":
ScipyTest().run()

```