1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE neupp.h.
6 Interface to ARPACK subroutines dneupd and sneupd.
7
8 ARPACK Authors
9 Richard Lehoucq
10 Danny Sorensen
11 Chao Yang
12 Dept. of Computational & Applied Mathematics
13 Rice University
14 Houston, Texas
15 */
16
17 #ifndef NEUPP_H
18 #define NEUPP_H
19
20 #include <cstddef>
21 #include <string>
22 #include "arch.h"
23 #include "arpackf.h"
24
neupp(bool rvec,char HowMny,double dr[],double di[],double Z[],ARint ldz,double sigmar,double sigmai,double workv[],char bmat,ARint n,const std::string & which,ARint nev,double tol,double resid[],ARint ncv,double V[],ARint ldv,ARint iparam[],ARint ipntr[],double workd[],double workl[],ARint lworkl,ARint & info)25 inline void neupp(bool rvec, char HowMny, double dr[],
26 double di[], double Z[], ARint ldz, double sigmar,
27 double sigmai, double workv[], char bmat, ARint n,
28 const std::string& which, ARint nev, double tol, double resid[],
29 ARint ncv, double V[], ARint ldv, ARint iparam[],
30 ARint ipntr[], double workd[], double workl[],
31 ARint lworkl, ARint& info)
32
33 /*
34 c++ version of ARPACK routine dneupd.
35 This subroutine returns the converged approximations to eigenvalues
36 of A*z = lambda*B*z and (optionally):
37
38 (1) the corresponding approximate eigenvectors,
39 (2) an orthonormal basis for the associated approximate
40 invariant subspace,
41
42 There is negligible additional cost to obtain eigenvectors. An
43 orthonormal basis is always computed. There is an additional storage cost
44 of n*nev if both are requested (in this case a separate array Z must be
45 supplied).
46 The approximate eigenvalues and eigenvectors of A*z = lambda*B*z
47 are derived from approximate eigenvalues and eigenvectors of
48 of the linear operator OP prescribed by the MODE selection in the
49 call to naupp. naupp must be called before this routine is called.
50 These approximate eigenvalues and vectors are commonly called Ritz
51 values and Ritz vectors respectively. They are referred to as such
52 in the comments that follow. The computed orthonormal basis for the
53 invariant subspace corresponding to these Ritz values is referred to
54 as a Schur basis.
55 See documentation in the header of the subroutine naupp for
56 definition of OP as well as other terms and the relation of computed
57 Ritz values and Ritz vectors of OP with respect to the given problem
58 A*z = lambda*B*z. For a brief description, see definitions of
59 iparam[7], MODE and which in the documentation of naupp.
60
61 Parameters:
62
63 rvec (Input) Specifies whether Ritz vectors corresponding to the
64 Ritz value approximations to the eigenproblem A*z = lambda*B*z
65 are computed.
66 rvec = false: Compute Ritz values only.
67 rvec = true : Compute the Ritz vectors or Schur vectors.
68 See Remarks below.
69 HowMny (Input) Specifies the form of the basis for the invariant
70 subspace corresponding to the converged Ritz values that
71 is to be computed.
72 = 'A': Compute nev Ritz vectors;
73 = 'P': Compute nev Schur vectors;
74 dr (Output) Array of dimension nev+1.
75 If iparam[7] = 1,2 or 3 and sigmai=0.0 then on exit: dr
76 contains the real part of the Ritz approximations to the
77 eigenvalues of A*z = lambda*B*z.
78 If iparam[7] = 3, 4 and sigmai is not equal to zero, then on
79 exit: dr contains the real part of the Ritz values of OP
80 computed by naupp. A further computation must be performed by
81 the user to transform the Ritz values computed for OP by naupp
82 to those of the original system A*z = lambda*B*z. See remark 3.
83 di (Output) Array of dimension nev+1.
84 On exit, di contains the imaginary part of the Ritz value
85 approximations to the eigenvalues of A*z = lambda*B*z
86 associated with dr.
87 NOTE: When Ritz values are complex, they will come in complex
88 conjugate pairs. If eigenvectors are requested, the
89 corresponding Ritz vectors will also come in conjugate
90 pairs and the real and imaginary parts of these are
91 represented in two consecutive columns of the array Z
92 (see below).
93 Z (Output) Array of dimension nev*n if rvec = TRUE and HowMny =
94 'A'. if rvec = TRUE. and HowMny = 'A', then the contains
95 approximate eigenvectors (Ritz vectors) corresponding to the
96 NCONV=iparam[5] Ritz values for eigensystem A*z = lambda*B*z.
97 The complex Ritz vector associated with the Ritz value
98 with positive imaginary part is stored in two consecutive
99 columns. The first column holds the real part of the Ritz
100 vector and the second column holds the imaginary part. The
101 Ritz vector associated with the Ritz value with negative
102 imaginary part is simply the complex conjugate of the Ritz
103 vector associated with the positive imaginary part.
104 If rvec = .FALSE. or HowMny = 'P', then Z is not referenced.
105 NOTE: If if rvec = .TRUE. and a Schur basis is not required,
106 the array Z may be set equal to first nev+1 columns of
107 the Arnoldi basis array V computed by naupp. In this
108 case the Arnoldi basis will be destroyed and overwritten
109 with the eigenvector basis.
110 ldz (Input) Dimension of the vectors contained in Z. This
111 parameter MUST be set to n.
112 sigmar (Input) If iparam[7] = 3 or 4, represents the real part of
113 the shift. Not referenced if iparam[7] = 1 or 2.
114 sigmai (Input) If iparam[7] = 3 or 4, represents the imaginary part
115 of the shift. Not referenced if iparam[7] = 1 or 2. See
116 remark 3 below.
117 workv (Workspace) Array of dimension 3*ncv.
118 V (Input/Output) Array of dimension n*ncv+1.
119 Upon Input: V contains the ncv vectors of the Arnoldi basis
120 for OP as constructed by naupp.
121 Upon Output: If rvec = TRUE the first NCONV=iparam[5] columns
122 contain approximate Schur vectors that span the
123 desired invariant subspace. See Remark 2 below.
124 NOTE: If the array Z has been set equal to first nev+1 columns
125 of the array V and rvec = TRUE. and HowMny = 'A', then
126 the Arnoldi basis held by V has been overwritten by the
127 desired Ritz vectors. If a separate array Z has been
128 passed then the first NCONV=iparam[5] columns of V will
129 contain approximate Schur vectors that span the desired
130 invariant subspace.
131 workl (Input / Output) Array of length lworkl+1.
132 workl[1:ncv*ncv+3*ncv] contains information obtained in
133 naupp. They are not changed by neupp.
134 workl[ncv*ncv+3*ncv+1:3*ncv*ncv+6*ncv] holds the real and
135 imaginary part of the untransformed Ritz values, the upper
136 quasi-triangular matrix for H, and the associated matrix
137 representation of the invariant subspace for H.
138 ipntr (Input / Output) Array of length 14. Pointer to mark the
139 starting locations in the workl array for matrices/vectors
140 used by naupp and neupp.
141 ipntr[9]: pointer to the real part of the ncv RITZ values
142 of the original system.
143 ipntr[10]: pointer to the imaginary part of the ncv RITZ
144 values of the original system.
145 ipntr[11]: pointer to the ncv corresponding error bounds.
146 ipntr[12]: pointer to the ncv by ncv upper quasi-triangular
147 Schur matrix for H.
148 ipntr[13]: pointer to the ncv by ncv matrix of eigenvectors
149 of the upper Hessenberg matrix H. Only referenced
150 by neupp if rvec = TRUE. See Remark 2 below.
151 info (Output) Error flag.
152 = 0 : Normal exit.
153 = 1 : The Schur form computed by LAPACK routine dlahqr
154 could not be reordered by LAPACK routine dtrsen.
155 Re-enter subroutine neupp with iparam[5] = ncv and
156 increase the size of the arrays DR and DI to have
157 dimension at least dimension ncv and allocate at least
158 ncv columns for Z. NOTE: Not necessary if Z and V share
159 the same space. Please notify the authors if this error
160 occurs.
161 = -1 : n must be positive.
162 = -2 : nev must be positive.
163 = -3 : ncv must satisfy nev+2 <= ncv <= n.
164 = -5 : which must be one of 'LM','SM','LR','SR','LI','SI'.
165 = -6 : bmat must be one of 'I' or 'G'.
166 = -7 : Length of private work workl array is not sufficient.
167 = -8 : Error return from calculation of a real Schur form.
168 Informational error from LAPACK routine dlahqr.
169 = -9 : Error return from calculation of eigenvectors.
170 Informational error from LAPACK routine dtrevc.
171 = -10: iparam[7] must be 1,2,3,4.
172 = -11: iparam[7] = 1 and bmat = 'G' are incompatible.
173 = -12: HowMny = 'S' not yet implemented
174 = -13: HowMny must be one of 'A' or 'P' if rvec = TRUE.
175 = -14: naupp did not find any eigenvalues to sufficient
176 accuracy.
177
178 NOTE: The following arguments
179
180 bmat, n, which, nev, tol, resid, ncv, V, ldv, iparam,
181 ipntr, workd, workl, lworkl, info
182
183 must be passed directly to neupp following the last call
184 to naupp. These arguments MUST NOT BE MODIFIED between
185 the the last call to naupp and the call to neupp.
186
187 Remarks
188 1. Currently only HowMny = 'A' and 'P' are implemented.
189 2. Schur vectors are an orthogonal representation for the basis of
190 Ritz vectors. Thus, their numerical properties are often superior.
191 Let X' denote the transpose of X. If rvec = .TRUE. then the
192 relationship A * V[:,1:iparam[5]] = V[:,1:iparam[5]] * T, and
193 V[:,1:iparam[5]]' * V[:,1:iparam[5]] = I are approximately satisfied.
194 Here T is the leading submatrix of order iparam[5] of the real
195 upper quasi-triangular matrix stored workl[ipntr[12]]. That is,
196 T is block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
197 each 2-by-2 diagonal block has its diagonal elements equal and its
198 off-diagonal elements of opposite sign. Corresponding to each
199 2-by-2 diagonal block is a complex conjugate pair of Ritz values.
200 The real Ritz values are stored on the diagonal of T.
201 3. If iparam[7] = 3 or 4 and sigmai is not equal zero, then the user
202 must form the iparam[5] Rayleigh quotients in order to transform the
203 Ritz values computed by naupp for OP to those of A*z = lambda*B*z.
204 Set rvec = TRUE. and HowMny = 'A', and compute
205 Z[:,I]' * A * Z[:,I] if di[I] = 0.
206 If di[I] is not equal to zero and di[I+1] = - D[I],
207 then the desired real and imaginary parts of the Ritz value are
208 Z[:,I]' * A * Z[:,I] + Z[:,I+1]' * A * Z[:,I+1],
209 Z[:,I]' * A * Z[:,I+1] - Z[:,I+1]' * A * Z[:,I], respectively.
210 Another possibility is to set rvec = .true. and HowMny = 'P' and
211 compute V[:,1:iparam[5]]' * A * V[:,1:iparam[5]] and then an upper
212 quasi-triangular matrix of order iparam[5] is computed. See remark
213 2 above.
214 */
215
216 {
217
218 ARint irvec;
219 ARlogical* iselect;
220 double* iZ;
221
222 irvec = (ARint) rvec;
223 iselect = new ARlogical[ncv];
224 iZ = (Z == NULL) ? &V[1] : Z;
225
226 F77NAME(dneupd)(&irvec, &HowMny, iselect, dr, di, iZ, &ldz, &sigmar,
227 &sigmai, &workv[1], &bmat, &n, which.c_str(), &nev, &tol,
228 resid, &ncv, &V[1], &ldv, &iparam[1], &ipntr[1],
229 &workd[1], &workl[1], &lworkl, &info);
230
231 delete[] iselect;
232
233 } // neupp (double).
234
neupp(bool rvec,char HowMny,float dr[],float di[],float Z[],ARint ldz,float sigmar,float sigmai,float workv[],char bmat,ARint n,const std::string & which,ARint nev,float tol,float resid[],ARint ncv,float V[],ARint ldv,ARint iparam[],ARint ipntr[],float workd[],float workl[],ARint lworkl,ARint & info)235 inline void neupp(bool rvec, char HowMny, float dr[],
236 float di[], float Z[], ARint ldz, float sigmar,
237 float sigmai, float workv[], char bmat, ARint n,
238 const std::string& which, ARint nev, float tol, float resid[],
239 ARint ncv, float V[], ARint ldv, ARint iparam[],
240 ARint ipntr[], float workd[], float workl[],
241 ARint lworkl, ARint& info)
242
243 /*
244 c++ version of ARPACK routine sneupd. The only difference between
245 sneupd and dneupd is that in the former function all vectors have
246 single precision elements and in the latter all vectors have double
247 precision elements.
248 */
249
250 {
251
252 ARint irvec;
253 ARlogical* iselect;
254 float* iZ;
255
256 irvec = (ARint) rvec;
257 iselect = new ARlogical[ncv];
258 iZ = (Z == NULL) ? &V[1] : Z;
259
260 F77NAME(sneupd)(&irvec, &HowMny, iselect, dr, di, iZ, &ldz, &sigmar,
261 &sigmai, &workv[1], &bmat, &n, which.c_str(), &nev, &tol,
262 resid, &ncv, &V[1], &ldv, &iparam[1], &ipntr[1],
263 &workd[1], &workl[1], &lworkl, &info );
264
265 delete[] iselect;
266
267 } // neupp (float).
268
269 #endif // NEUPP_H
270
271