1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ceupp.h.
6 Interface to ARPACK subroutines zneupd and cneupd.
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 CEUPP_H
18 #define CEUPP_H
19
20 #include <cstddef>
21 #include <string>
22 #include "arch.h"
23 #include "arpackf.h"
24
ceupp(bool rvec,char HowMny,arcomplex<double> d[],arcomplex<double> Z[],ARint ldz,arcomplex<double> sigma,arcomplex<double> workev[],char bmat,ARint n,const std::string & which,ARint nev,double tol,arcomplex<double> resid[],ARint ncv,arcomplex<double> V[],ARint ldv,ARint iparam[],ARint ipntr[],arcomplex<double> workd[],arcomplex<double> workl[],ARint lworkl,double rwork[],ARint & info)25 inline void ceupp(bool rvec, char HowMny, arcomplex<double> d[],
26 arcomplex<double> Z[], ARint ldz, arcomplex<double> sigma,
27 arcomplex<double> workev[], char bmat, ARint n, const std::string& which,
28 ARint nev, double tol, arcomplex<double> resid[], ARint ncv,
29 arcomplex<double> V[], ARint ldv, ARint iparam[],
30 ARint ipntr[], arcomplex<double> workd[],
31 arcomplex<double> workl[], ARint lworkl, double rwork[],
32 ARint& info)
33
34 /*
35 c++ version of ARPACK routine zneupd.
36 This subroutine returns the converged approximations to eigenvalues
37 of A*z = lambda*B*z and (optionally):
38
39 (1) the corresponding approximate eigenvectors,
40 (2) an orthonormal basis for the associated approximate
41 invariant subspace,
42
43 There is negligible additional cost to obtain eigenvectors. An
44 orthonormal basis is always computed. There is an additional storage cost
45 of n*nev if both are requested (in this case a separate array Z must be
46 supplied).
47 The approximate eigenvalues and eigenvectors of A*z = lambda*B*z
48 are derived from approximate eigenvalues and eigenvectors of
49 of the linear operator OP prescribed by the MODE selection in the
50 call to caupp. caupp must be called before this routine is called.
51 These approximate eigenvalues and vectors are commonly called Ritz
52 values and Ritz vectors respectively. They are referred to as such
53 in the comments that follow. The computed orthonormal basis for the
54 invariant subspace corresponding to these Ritz values is referred to
55 as a Schur basis.
56 See documentation in the header of the subroutine caupp for
57 definition of OP as well as other terms and the relation of computed
58 Ritz values and Ritz vectors of OP with respect to the given problem
59 A*z = lambda*B*z. For a brief description, see definitions of
60 iparam[7], MODE and which in the documentation of caupp.
61
62 Parameters:
63
64 rvec (Input) Specifies whether a basis for the invariant subspace
65 corresponding to the converged Ritz value approximations for
66 the eigenproblem A*z = lambda*B*z is computed.
67 rvec = false: Compute Ritz values only.
68 rvec = true : Compute the Ritz vectors or Schur vectors.
69 See Remarks below.
70 HowMny (Input) Specifies the form of the basis for the invariant
71 subspace corresponding to the converged Ritz values that
72 is to be computed.
73 = 'A': Compute nev Ritz vectors;
74 = 'P': Compute nev Schur vectors;
75 d (Output) Array of dimension nev+1. D contains the Ritz
76 approximations to the eigenvalues lambda for A*z = lambda*B*z.
77 Z (Output) Array of dimension nev*n. If rvec = TRUE. and
78 HowMny = 'A', then Z contains approximate eigenvectors (Ritz
79 vectors) corresponding to the NCONV=iparam[5] Ritz values for
80 eigensystem A*z = lambda*B*z.
81 If rvec = .FALSE. or HowMny = 'P', then Z is not referenced.
82 NOTE: If if rvec = .TRUE. and a Schur basis is not required,
83 the array Z may be set equal to first nev+1 columns of
84 the Arnoldi basis array V computed by caupp. In this
85 case the Arnoldi basis will be destroyed and overwritten
86 with the eigenvector basis.
87 ldz (Input) Dimension of the vectors contained in Z. This
88 parameter MUST be set to n.
89 sigma (Input) If iparam[7] = 3, sigma represents the shift. Not
90 referenced if iparam[7] = 1 or 2.
91 workv (Workspace) Array of dimension 2*ncv.
92 V (Input/Output) Array of dimension n*ncv+1.
93 Upon Input: V contains the ncv vectors of the Arnoldi basis
94 for OP as constructed by caupp.
95 Upon Output: If rvec = TRUE the first NCONV=iparam[5] columns
96 contain approximate Schur vectors that span the
97 desired invariant subspace.
98 NOTE: If the array Z has been set equal to first nev+1 columns
99 of the array V and rvec = TRUE. and HowMny = 'A', then
100 the Arnoldi basis held by V has been overwritten by the
101 desired Ritz vectors. If a separate array Z has been
102 passed then the first NCONV=iparam[5] columns of V will
103 contain approximate Schur vectors that span the desired
104 invariant subspace.
105 workl (Input / Output) Array of length lworkl+1.
106 workl[1:ncv*ncv+3*ncv] contains information obtained in
107 caupp. They are not changed by ceupp.
108 workl[ncv*ncv+3*ncv+1:3*ncv*ncv+4*ncv] holds the untransformed
109 Ritz values, the untransformed error estimates of the Ritz
110 values, the upper triangular matrix for H, and the associated
111 matrix representation of the invariant subspace for H.
112 ipntr (Input / Output) Array of length 14. Pointer to mark the
113 starting locations in the workl array for matrices/vectors
114 used by caupp and ceupp.
115 ipntr[9]: pointer to the ncv RITZ values of the original
116 system.
117 ipntr[11]: pointer to the ncv corresponding error estimates.
118 ipntr[12]: pointer to the ncv by ncv upper triangular
119 Schur matrix for H.
120 ipntr[13]: pointer to the ncv by ncv matrix of eigenvectors
121 of the upper Hessenberg matrix H. Only referenced
122 by ceupp if rvec = TRUE. See Remark 2 below.
123 info (Output) Error flag.
124 = 0 : Normal exit.
125 = 1 : The Schur form computed by LAPACK routine csheqr
126 could not be reordered by LAPACK routine ztrsen.
127 Re-enter subroutine ceupp with iparam[5] = ncv and
128 increase the size of the array D to have
129 dimension at least dimension ncv and allocate at least
130 ncv columns for Z. NOTE: Not necessary if Z and V share
131 the same space. Please notify the authors if this error
132 occurs.
133 = -1 : n must be positive.
134 = -2 : nev must be positive.
135 = -3 : ncv must satisfy nev+1 <= ncv <= n.
136 = -5 : which must be one of 'LM','SM','LR','SR','LI','SI'.
137 = -6 : bmat must be one of 'I' or 'G'.
138 = -7 : Length of private work workl array is not sufficient.
139 = -8 : Error return from LAPACK eigenvalue calculation.
140 This should never happened.
141 = -9 : Error return from calculation of eigenvectors.
142 Informational error from LAPACK routine ztrevc.
143 = -10: iparam[7] must be 1, 2 or 3.
144 = -11: iparam[7] = 1 and bmat = 'G' are incompatible.
145 = -12: HowMny = 'S' not yet implemented.
146 = -13: HowMny must be one of 'A' or 'P' if rvec = TRUE.
147 = -14: caupp did not find any eigenvalues to sufficient
148 accuracy.
149
150 NOTE: The following arguments
151
152 bmat, n, which, nev, tol, resid, ncv, V, ldv, iparam,
153 ipntr, workd, workl, lworkl, rwork, info
154
155 must be passed directly to ceupp following the last call
156 to caupp. These arguments MUST NOT BE MODIFIED between
157 the the last call to caupp and the call to ceupp.
158
159 Remarks
160 1. Currently only HowMny = 'A' and 'P' are implemented.
161 2. Schur vectors are an orthogonal representation for the basis of
162 Ritz vectors. Thus, their numerical properties are often superior.
163 Let X' denote the transpose of X. If rvec = .TRUE. then the
164 relationship A * V[:,1:iparam[5]] = V[:,1:iparam[5]] * T, and
165 V[:,1:iparam[5]]' * V[:,1:iparam[5]] = I are approximately satisfied.
166 Here T is the leading submatrix of order iparam[5] of the real
167 upper quasi-triangular matrix stored workl[ipntr[12]].
168 */
169
170 {
171
172 ARint irvec;
173 ARlogical* iselect;
174 arcomplex<double>* iZ;
175
176 irvec = (ARint) rvec;
177 iselect = new ARlogical[ncv];
178 iZ = (Z == NULL) ? &V[1] : Z;
179
180 F77NAME(zneupd)(&irvec, &HowMny, iselect, d, iZ, &ldz, &sigma,
181 &workev[1], &bmat, &n, which.c_str(), &nev, &tol, resid,
182 &ncv, &V[1], &ldv, &iparam[1], &ipntr[1],
183 &workd[1], &workl[1], &lworkl, &rwork[1], &info);
184
185 delete[] iselect;
186
187 } // ceupp (arcomplex<double>).
188
ceupp(bool rvec,char HowMny,arcomplex<float> d[],arcomplex<float> Z[],ARint ldz,arcomplex<float> sigma,arcomplex<float> workev[],char bmat,ARint n,const std::string & which,ARint nev,float tol,arcomplex<float> resid[],ARint ncv,arcomplex<float> V[],ARint ldv,ARint iparam[],ARint ipntr[],arcomplex<float> workd[],arcomplex<float> workl[],ARint lworkl,float rwork[],ARint & info)189 inline void ceupp(bool rvec, char HowMny, arcomplex<float> d[],
190 arcomplex<float> Z[], ARint ldz, arcomplex<float> sigma,
191 arcomplex<float> workev[], char bmat, ARint n, const std::string& which,
192 ARint nev, float tol, arcomplex<float> resid[], ARint ncv,
193 arcomplex<float> V[], ARint ldv, ARint iparam[],
194 ARint ipntr[], arcomplex<float> workd[],
195 arcomplex<float> workl[], ARint lworkl, float rwork[],
196 ARint& info)
197
198 /*
199 c++ version of ARPACK routine cneupd. The only difference between
200 cneupd and zneupd is that in the former function all vectors have
201 single precision elements and in the latter all vectors have double
202 precision elements.
203 */
204
205 {
206
207 ARint irvec;
208 ARlogical* iselect;
209 arcomplex<float>* iZ;
210
211 irvec = (ARint) rvec;
212 iselect = new ARlogical[ncv];
213 iZ = (Z == NULL) ? &V[1] : Z;
214
215 F77NAME(cneupd)(&irvec, &HowMny, iselect, d, iZ, &ldz, &sigma,
216 &workev[1], &bmat, &n, which.c_str(), &nev, &tol, resid,
217 &ncv, &V[1], &ldv, &iparam[1], &ipntr[1],
218 &workd[1], &workl[1], &lworkl, &rwork[1], &info);
219
220 delete[] iselect;
221
222 } // ceupp (arcomplex<float>).
223
224 #endif // CEUPP_H
225