1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
|
import numpy as N
import _arpack
import warnings
__all___=['ArpackException','ARPACK_eigs', 'ARPACK_gen_eigs']
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+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
# Storage for the Arnoldi basis vectors
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
_arpack.debug.mnaitr = 0
_arpack.debug.mnapps = 0
_arpack.debug.mnaupd = 1
_arpack.debug.mnaup2 = 0
_arpack.debug.mneigh = 0
_arpack.debug.mneupd = 1
def init_postproc_workspace(n, nev, ncv):
# Used as workspace and to return eigenvectors if requested. Not touched if
# eigenvectors are not requested
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, 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 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
eigenvalues, nev eigenvalues close to sigma are calculated. The user needs
to provide routines that calculate [M]*x and solve [A]-sigma*[M]*x = b for x.
Arguments
=========
matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> [M]*x
sigma_solve -- sigma_solve(b) -> x, where [A]-sigma*[M]*x = b
n -- Matrix dimension of the problem
sigma -- Eigenvalue spectral shift real value
nev -- Number of eigenvalues to calculate
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
=============
(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.
"""
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)
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
# Some random arpack parameter (I think it tells ARPACK to solve the
# 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
# 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 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)
elif ido == 1: # Perform y = inv[A - sigma*M]*M*x using saved M*x
# Mx = [M]*x where it was saved by ARPACK
Mx = workd[ipntr[2]-1:ipntr[2]+n-1]
workd[ipntr[1]-1:ipntr[1]+n-1] = sigma_solve(Mx)
elif ido == 2: # Perform y = M*x
x = workd[ipntr[0]-1:ipntr[0]+n-1]
workd[ipntr[1]-1:ipntr[1]+n-1] = matvec(x)
else: # Finished, or error
break
if info == 1:
warn.warn("Maximum number of iterations taken: %s"%iparam[2])
elif info != 0:
raise ArpackException(info)
return (ncv, resid, iparam, ipntr, v, workd, workl, info)
|
