4
__all___=['ArpackException','ARPACK_eigs', 'ARPACK_gen_eigs']
6
class ArpackException(RuntimeError):
7
ARPACKErrors = { 0: """Normal exit.""",
8
3: """No shifts could be applied during a cycle of the
9
Implicitly restarted Arnoldi iteration. One possibility
10
is to increase the size of NCV relative to NEV.""",
11
-1: """N must be positive.""",
12
-2: """NEV must be positive.""",
13
-3: """NCV-NEV >= 2 and less than or equal to N.""",
14
-4: """The maximum number of Arnoldi update iteration
15
must be greater than zero.""",
16
-5: """WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'""",
17
-6: """BMAT must be one of 'I' or 'G'.""",
18
-7: """Length of private work array is not sufficient.""",
19
-8: """Error return from LAPACK eigenvalue calculation;""",
20
-9: """Starting vector is zero.""",
21
-10: """IPARAM(7) must be 1,2,3,4.""",
22
-11: """IPARAM(7) = 1 and BMAT = 'G' are incompatable.""",
23
-12: """IPARAM(1) must be equal to 0 or 1.""",
24
-9999: """Could not build an Arnoldi factorization.
25
IPARAM(5) returns the size of the current Arnoldi
28
def __init__(self, info):
31
try: return self.ARPACKErrors[self.info]
32
except KeyError: return "Unknown ARPACK error"
34
def check_init(n, nev, ncv):
35
assert(nev <= n-4) # ARPACK seems to cause a segfault otherwise
37
ncv = min(2*nev+1, n-1)
38
maxitr = max(n, 1000) # Maximum number of iterations
41
def init_workspaces(n,nev,ncv):
42
ipntr = np.zeros(14, np.int32) # Pointers into memory structure used by F77 calls
43
d = np.zeros((ncv, 3), np.float64, order='FORTRAN') # Temp workspace
44
# Temp workspace/error residuals upon iteration completion
45
resid = np.zeros(n, np.float64)
46
workd = np.zeros(3*n, np.float64) # workspace
47
workl = np.zeros(3*ncv*ncv+6*ncv, np.float64) # more workspace
48
# Storage for the Arnoldi basis vectors
49
v = np.zeros((n, ncv), dtype=np.float64, order='FORTRAN')
50
return (ipntr, d, resid, workd, workl, v)
53
# Causes various debug info to be printed by ARPACK
54
_arpack.debug.ndigit = -3
55
_arpack.debug.logfil = 6
56
_arpack.debug.mnaitr = 0
57
_arpack.debug.mnapps = 0
58
_arpack.debug.mnaupd = 1
59
_arpack.debug.mnaup2 = 0
60
_arpack.debug.mneigh = 0
61
_arpack.debug.mneupd = 1
63
def init_postproc_workspace(n, nev, ncv):
64
# Used as workspace and to return eigenvectors if requested. Not touched if
65
# eigenvectors are not requested
66
workev = np.zeros(3*ncv, np.float64) # More workspace
67
select = np.zeros(ncv, np.int32) # Used as internal workspace since dneupd
68
# parameter HOWMNY == 'A'
69
return (workev, select)
71
def postproc(n, nev, ncv, sigmar, sigmai, bmat, which,
72
tol, resid, v, iparam, ipntr, workd, workl, info):
73
workev, select = init_postproc_workspace(n, nev, ncv)
75
# Postprocess the Arnouldi vectors to extract eigenvalues/vectors
76
# If dneupd's first paramter is 'True' the eigenvectors are also calculated,
77
# 'False' only the eigenvalues
78
dr,di,z,info = _arpack.dneupd(
79
True, 'A', select, sigmar, sigmai, workev, bmat, which, nev, tol, resid, v,
80
iparam, ipntr, workd, workl, info)
82
if np.abs(di[:-1]).max() == 0: dr = dr[:-1]
83
else: dr = dr[:-1] + 1j*di[:-1]
87
def ARPACK_eigs(matvec, n, nev, which='SM', ncv=None, tol=1e-14):
89
Calculate eigenvalues for system with matrix-vector product matvec, dimension n
93
matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> A*x
94
n -- Matrix dimension of the problem
95
nev -- Number of eigenvalues to calculate
96
which -- Spectrum selection. See details below. Defaults to 'SM'
97
ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev+1
98
tol -- Numerical tollerance for Arnouldi iteration convergence. Defaults to 1e-14
102
which can take one of several values:
104
'LM' -> Request eigenvalues with largest magnitude.
105
'SM' -> Request eigenvalues with smallest magnitude.
106
'LR' -> Request eigenvalues with largest real part.
107
'SR' -> Request eigenvalues with smallest real part.
108
'LI' -> Request eigenvalues with largest imaginary part.
109
'SI' -> Request eigenvalues with smallest imaginary part.
113
(eig_vals, eig_vecs) where eig_vals are the requested eigenvalues and
114
eig_vecs the corresponding eigenvectors. If all the eigenvalues are real,
115
eig_vals is a real array but if some eigenvalues are complex it is a
119
bmat = 'I' # Standard eigenproblem
120
ncv, resid, iparam, ipntr, v, workd, workl, info = ARPACK_iteration(
121
matvec, lambda x: x, n, bmat, which, nev, tol, ncv, mode=1)
122
return postproc(n, nev, ncv, 0., 0., bmat, which, tol,
123
resid, v, iparam, ipntr, workd, workl, info)
125
def ARPACK_gen_eigs(matvec, sigma_solve, n, sigma, nev, which='LR', ncv=None, tol=1e-14):
127
Calculate eigenvalues close to sigma for generalised eigen system
129
Given a system [A]x = k_i*[M]x where [A] and [M] are matrices and k_i are
130
eigenvalues, nev eigenvalues close to sigma are calculated. The user needs
131
to provide routines that calculate [M]*x and solve [A]-sigma*[M]*x = b for x.
135
matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> [M]*x
136
sigma_solve -- sigma_solve(b) -> x, where [A]-sigma*[M]*x = b
137
n -- Matrix dimension of the problem
138
sigma -- Eigenvalue spectral shift real value
139
nev -- Number of eigenvalues to calculate
140
which -- Spectrum selection. See details below. Defaults to 'LR'
141
ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev+1
142
tol -- Numerical tollerance for Arnouldi iteration convergence. Defaults to 1e-14
147
The spectrum of the orignal system is shifted by sigma. This transforms the
148
original eigenvalues to be 1/(original_eig-sigma) in the shifted
149
system. ARPACK then operates on the shifted system, transforming it back to
150
the original system in a postprocessing step.
152
The spectrum shift causes eigenvalues close to sigma to become very large
153
in the transformed system. This allows quick convergence for these
154
eigenvalues. This is particularly useful if a system has a number of
155
trivial zero-eigenvalues that are to be ignored.
159
which can take one of several values:
161
'LM' -> Request spectrum shifted eigenvalues with largest magnitude.
162
'SM' -> Request spectrum shifted eigenvalues with smallest magnitude.
163
'LR' -> Request spectrum shifted eigenvalues with largest real part.
164
'SR' -> Request spectrum shifted eigenvalues with smallest real part.
165
'LI' -> Request spectrum shifted eigenvalues with largest imaginary part.
166
'SI' -> Request spectrum shifted eigenvalues with smallest imaginary part.
168
The effect on the actual system is:
169
'LM' -> Eigenvalues closest to sigma on the complex plane
170
'LR' -> Eigenvalues with real part > sigma, provided they exist
175
(eig_vals, eig_vecs) where eig_vals are the requested eigenvalues and
176
eig_vecs the corresponding eigenvectors. If all the eigenvalues are real,
177
eig_vals is a real array but if some eigenvalues are complex it is a
178
complex array. The eigenvalues and vectors correspond to the original
179
system, not the shifted system. The shifted system is only used interally.
182
bmat = 'G' # Generalised eigenproblem
183
ncv, resid, iparam, ipntr, v, workd, workl, info = ARPACK_iteration(
184
matvec, sigma_solve, n, bmat, which, nev, tol, ncv, mode=3)
187
return postproc(n, nev, ncv, sigmar, sigmai, bmat, which, tol,
188
resid, v, iparam, ipntr, workd, workl, info)
190
def ARPACK_iteration(matvec, sigma_solve, n, bmat, which, nev, tol, ncv, mode):
191
ncv, maxitr = check_init(n, nev, ncv)
192
ipntr, d, resid, workd, workl, v = init_workspaces(n,nev,ncv)
194
ishfts = 1 # Some random arpack parameter
195
# Some random arpack parameter (I think it tells ARPACK to solve the
196
# general eigenproblem using shift-invert
197
iparam = np.zeros(11, np.int32) # Array with assorted extra paramters for F77 call
198
iparam[[0,2,6]] = ishfts, maxitr, mode
199
ido = 0 # Communication variable used by ARPACK to tell the user what to do
200
info = 0 # Used for error reporting
201
# Arnouldi iteration.
203
ido,resid,v,iparam,ipntr,info = _arpack.dnaupd(
204
ido, bmat, which, nev, tol, resid, v, iparam, ipntr, workd, workl, info)
205
if ido == -1 or ido == 1 and mode not in (3,4):
206
# Perform y = inv[A - sigma*M]*M*x
207
x = workd[ipntr[0]-1:ipntr[0]+n-1]
208
Mx = matvec(x) # Mx = [M]*x
209
workd[ipntr[1]-1:ipntr[1]+n-1] = sigma_solve(Mx)
210
elif ido == 1: # Perform y = inv[A - sigma*M]*M*x using saved M*x
211
# Mx = [M]*x where it was saved by ARPACK
212
Mx = workd[ipntr[2]-1:ipntr[2]+n-1]
213
workd[ipntr[1]-1:ipntr[1]+n-1] = sigma_solve(Mx)
214
elif ido == 2: # Perform y = M*x
215
x = workd[ipntr[0]-1:ipntr[0]+n-1]
216
workd[ipntr[1]-1:ipntr[1]+n-1] = matvec(x)
217
else: # Finished, or error
220
warn.warn("Maximum number of iterations taken: %s"%iparam[2])
222
raise ArpackException(info)
224
return (ncv, resid, iparam, ipntr, v, workd, workl, info)
added
removed
scipy/sparse/linalg/eigen/arpack/speigs.py
