[Scipy-svn] r2324 - in trunk/Lib/sandbox/arpack: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Nov 20 15:29:04 EST 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
-# A,B = pickle.load(file('/tmp/mats.pickle'))
+# A,B = pickle.load(file('mats.pickle'))
# 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()
More information about the Scipy-svn
mailing list