1 /* ./src_f77/dneupd.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Common Block Declarations */
9
10 struct {
11 integer logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, msapps,
12 msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets,
13 mneupd, mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd;
14 } debug_DNEUPD;
15
16 #define debug_1 debug_DNEUPD
17
18 struct {
19 integer nopx, nbx, nrorth, nitref, nrstrt;
20 real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd,
21 tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, tcaupd, tcaup2,
22 tcaitr, tceigh, tcgets, tcapps, tcconv, tmvopx, tmvbx, tgetv0,
23 titref, trvec;
24 } timing_DNEUPD;
25
26 #define timing_1 timing_DNEUPD
27
28 /* Table of constant values */
29
30 static doublereal c_b3 = .66666666666666663;
31 static integer c__1 = 1;
32 static doublereal c_b37 = 0.;
33 static doublereal c_b38 = 1.;
34 static logical c_true = TRUE_;
35 static doublereal c_b64 = -1.;
36
37 /* \BeginDoc */
38
39 /* \Name: dneupd */
40
41 /* \Description: */
42
43 /* This subroutine returns the converged approximations to eigenvalues */
44 /* of A*z = lambda*B*z and (optionally): */
45
46 /* (1) The corresponding approximate eigenvectors; */
47
48 /* (2) An orthonormal basis for the associated approximate */
49 /* invariant subspace; */
50
51 /* (3) Both. */
52
53 /* There is negligible additional cost to obtain eigenvectors. An orthonormal */
54 /* basis is always computed. There is an additional storage cost of n*nev */
55 /* if both are requested (in this case a separate array Z must be supplied). */
56
57 /* The approximate eigenvalues and eigenvectors of A*z = lambda*B*z */
58 /* are derived from approximate eigenvalues and eigenvectors of */
59 /* of the linear operator OP prescribed by the MODE selection in the */
60 /* call to DNAUPD . DNAUPD must be called before this routine is called. */
61 /* These approximate eigenvalues and vectors are commonly called Ritz */
62 /* values and Ritz vectors respectively. They are referred to as such */
63 /* in the comments that follow. The computed orthonormal basis for the */
64 /* invariant subspace corresponding to these Ritz values is referred to as a */
65 /* Schur basis. */
66
67 /* See documentation in the header of the subroutine DNAUPD for */
68 /* definition of OP as well as other terms and the relation of computed */
69 /* Ritz values and Ritz vectors of OP with respect to the given problem */
70 /* A*z = lambda*B*z. For a brief description, see definitions of */
71 /* IPARAM(7), MODE and WHICH in the documentation of DNAUPD . */
72
73 /* \Usage: */
74 /* call dneupd */
75 /* ( RVEC, HOWMNY, SELECT, DR, DI, Z, LDZ, SIGMAR, SIGMAI, WORKEV, BMAT, */
76 /* N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, WORKL, */
77 /* LWORKL, INFO ) */
78
79 /* \Arguments: */
80 /* RVEC LOGICAL (INPUT) */
81 /* Specifies whether a basis for the invariant subspace corresponding */
82 /* to the converged Ritz value approximations for the eigenproblem */
83 /* A*z = lambda*B*z is computed. */
84
85 /* RVEC = .FALSE. Compute Ritz values only. */
86
87 /* RVEC = .TRUE. Compute the Ritz vectors or Schur vectors. */
88 /* See Remarks below. */
89
90 /* HOWMNY Character*1 (INPUT) */
91 /* Specifies the form of the basis for the invariant subspace */
92 /* corresponding to the converged Ritz values that is to be computed. */
93
94 /* = 'A': Compute NEV Ritz vectors; */
95 /* = 'P': Compute NEV Schur vectors; */
96 /* = 'S': compute some of the Ritz vectors, specified */
97 /* by the logical array SELECT. */
98
99 /* SELECT Logical array of dimension NCV. (INPUT) */
100 /* If HOWMNY = 'S', SELECT specifies the Ritz vectors to be */
101 /* computed. To select the Ritz vector corresponding to a */
102 /* Ritz value (DR(j), DI(j)), SELECT(j) must be set to .TRUE.. */
103 /* If HOWMNY = 'A' or 'P', SELECT is used as internal workspace. */
104
105 /* DR Double precision array of dimension NEV+1. (OUTPUT) */
106 /* If IPARAM(7) = 1,2 or 3 and SIGMAI=0.0 then on exit: DR contains */
107 /* the real part of the Ritz approximations to the eigenvalues of */
108 /* A*z = lambda*B*z. */
109 /* If IPARAM(7) = 3, 4 and SIGMAI is not equal to zero, then on exit: */
110 /* DR contains the real part of the Ritz values of OP computed by */
111 /* DNAUPD . A further computation must be performed by the user */
112 /* to transform the Ritz values computed for OP by DNAUPD to those */
113 /* of the original system A*z = lambda*B*z. See remark 3 below. */
114
115 /* DI Double precision array of dimension NEV+1. (OUTPUT) */
116 /* On exit, DI contains the imaginary part of the Ritz value */
117 /* approximations to the eigenvalues of A*z = lambda*B*z associated */
118 /* with DR. */
119
120 /* NOTE: When Ritz values are complex, they will come in complex */
121 /* conjugate pairs. If eigenvectors are requested, the */
122 /* corresponding Ritz vectors will also come in conjugate */
123 /* pairs and the real and imaginary parts of these are */
124 /* represented in two consecutive columns of the array Z */
125 /* (see below). */
126
127 /* Z Double precision N by NEV+1 array if RVEC = .TRUE. and HOWMNY = 'A'. (OUTPUT) */
128 /* On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of */
129 /* Z represent approximate eigenvectors (Ritz vectors) corresponding */
130 /* to the NCONV=IPARAM(5) Ritz values for eigensystem */
131 /* A*z = lambda*B*z. */
132
133 /* The complex Ritz vector associated with the Ritz value */
134 /* with positive imaginary part is stored in two consecutive */
135 /* columns. The first column holds the real part of the Ritz */
136 /* vector and the second column holds the imaginary part. The */
137 /* Ritz vector associated with the Ritz value with negative */
138 /* imaginary part is simply the complex conjugate of the Ritz vector */
139 /* associated with the positive imaginary part. */
140
141 /* If RVEC = .FALSE. or HOWMNY = 'P', then Z is not referenced. */
142
143 /* NOTE: If if RVEC = .TRUE. and a Schur basis is not required, */
144 /* the array Z may be set equal to first NEV+1 columns of the Arnoldi */
145 /* basis array V computed by DNAUPD . In this case the Arnoldi basis */
146 /* will be destroyed and overwritten with the eigenvector basis. */
147
148 /* LDZ Integer. (INPUT) */
149 /* The leading dimension of the array Z. If Ritz vectors are */
150 /* desired, then LDZ >= max( 1, N ). In any case, LDZ >= 1. */
151
152 /* SIGMAR Double precision (INPUT) */
153 /* If IPARAM(7) = 3 or 4, represents the real part of the shift. */
154 /* Not referenced if IPARAM(7) = 1 or 2. */
155
156 /* SIGMAI Double precision (INPUT) */
157 /* If IPARAM(7) = 3 or 4, represents the imaginary part of the shift. */
158 /* Not referenced if IPARAM(7) = 1 or 2. See remark 3 below. */
159
160 /* WORKEV Double precision work array of dimension 3*NCV. (WORKSPACE) */
161
162 /* **** The remaining arguments MUST be the same as for the **** */
163 /* **** call to DNAUPD that was just completed. **** */
164
165 /* NOTE: The remaining arguments */
166
167 /* BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, */
168 /* WORKD, WORKL, LWORKL, INFO */
169
170 /* must be passed directly to DNEUPD following the last call */
171 /* to DNAUPD . These arguments MUST NOT BE MODIFIED between */
172 /* the the last call to DNAUPD and the call to DNEUPD . */
173
174 /* Three of these parameters (V, WORKL, INFO) are also output parameters: */
175
176 /* V Double precision N by NCV array. (INPUT/OUTPUT) */
177
178 /* Upon INPUT: the NCV columns of V contain the Arnoldi basis */
179 /* vectors for OP as constructed by DNAUPD . */
180
181 /* Upon OUTPUT: If RVEC = .TRUE. the first NCONV=IPARAM(5) columns */
182 /* contain approximate Schur vectors that span the */
183 /* desired invariant subspace. See Remark 2 below. */
184
185 /* NOTE: If the array Z has been set equal to first NEV+1 columns */
186 /* of the array V and RVEC=.TRUE. and HOWMNY= 'A', then the */
187 /* Arnoldi basis held by V has been overwritten by the desired */
188 /* Ritz vectors. If a separate array Z has been passed then */
189 /* the first NCONV=IPARAM(5) columns of V will contain approximate */
190 /* Schur vectors that span the desired invariant subspace. */
191
192 /* WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE) */
193 /* WORKL(1:ncv*ncv+3*ncv) contains information obtained in */
194 /* dnaupd . They are not changed by dneupd . */
195 /* WORKL(ncv*ncv+3*ncv+1:3*ncv*ncv+6*ncv) holds the */
196 /* real and imaginary part of the untransformed Ritz values, */
197 /* the upper quasi-triangular matrix for H, and the */
198 /* associated matrix representation of the invariant subspace for H. */
199
200 /* Note: IPNTR(9:13) contains the pointer into WORKL for addresses */
201 /* of the above information computed by dneupd . */
202 /* ------------------------------------------------------------- */
203 /* IPNTR(9): pointer to the real part of the NCV RITZ values of the */
204 /* original system. */
205 /* IPNTR(10): pointer to the imaginary part of the NCV RITZ values of */
206 /* the original system. */
207 /* IPNTR(11): pointer to the NCV corresponding error bounds. */
208 /* IPNTR(12): pointer to the NCV by NCV upper quasi-triangular */
209 /* Schur matrix for H. */
210 /* IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors */
211 /* of the upper Hessenberg matrix H. Only referenced by */
212 /* dneupd if RVEC = .TRUE. See Remark 2 below. */
213 /* ------------------------------------------------------------- */
214
215 /* INFO Integer. (OUTPUT) */
216 /* Error flag on output. */
217
218 /* = 0: Normal exit. */
219
220 /* = 1: The Schur form computed by LAPACK routine dlahqr */
221 /* could not be reordered by LAPACK routine dtrsen . */
222 /* Re-enter subroutine dneupd with IPARAM(5)=NCV and */
223 /* increase the size of the arrays DR and DI to have */
224 /* dimension at least dimension NCV and allocate at least NCV */
225 /* columns for Z. NOTE: Not necessary if Z and V share */
226 /* the same space. Please notify the authors if this error */
227 /* occurs. */
228
229 /* = -1: N must be positive. */
230 /* = -2: NEV must be positive. */
231 /* = -3: NCV-NEV >= 2 and less than or equal to N. */
232 /* = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI' */
233 /* = -6: BMAT must be one of 'I' or 'G'. */
234 /* = -7: Length of private work WORKL array is not sufficient. */
235 /* = -8: Error return from calculation of a real Schur form. */
236 /* Informational error from LAPACK routine dlahqr . */
237 /* = -9: Error return from calculation of eigenvectors. */
238 /* Informational error from LAPACK routine dtrevc . */
239 /* = -10: IPARAM(7) must be 1,2,3,4. */
240 /* = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible. */
241 /* = -12: HOWMNY = 'S' not yet implemented */
242 /* = -13: HOWMNY must be one of 'A' or 'P' if RVEC = .true. */
243 /* = -14: DNAUPD did not find any eigenvalues to sufficient */
244 /* accuracy. */
245 /* = -15: DNEUPD got a different count of the number of converged */
246 /* Ritz values than DNAUPD got. This indicates the user */
247 /* probably made an error in passing data from DNAUPD to */
248 /* DNEUPD or that the data was modified before entering */
249 /* DNEUPD */
250
251 /* \BeginLib */
252
253 /* \References: */
254 /* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
255 /* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
256 /* pp 357-385. */
257 /* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
258 /* Restarted Arnoldi Iteration", Rice University Technical Report */
259 /* TR95-13, Department of Computational and Applied Mathematics. */
260 /* 3. B.N. Parlett & Y. Saad, "Complex Shift and Invert Strategies for */
261 /* Real Matrices", Linear Algebra and its Applications, vol 88/89, */
262 /* pp 575-595, (1987). */
263
264 /* \Routines called: */
265 /* ivout ARPACK utility routine that prints integers. */
266 /* dmout ARPACK utility routine that prints matrices */
267 /* dvout ARPACK utility routine that prints vectors. */
268 /* dgeqr2 LAPACK routine that computes the QR factorization of */
269 /* a matrix. */
270 /* dlacpy LAPACK matrix copy routine. */
271 /* dlahqr LAPACK routine to compute the real Schur form of an */
272 /* upper Hessenberg matrix. */
273 /* dlamch LAPACK routine that determines machine constants. */
274 /* dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully. */
275 /* dlaset LAPACK matrix initialization routine. */
276 /* dorm2r LAPACK routine that applies an orthogonal matrix in */
277 /* factored form. */
278 /* dtrevc LAPACK routine to compute the eigenvectors of a matrix */
279 /* in upper quasi-triangular form. */
280 /* dtrsen LAPACK routine that re-orders the Schur form. */
281 /* dtrmm Level 3 BLAS matrix times an upper triangular matrix. */
282 /* dger Level 2 BLAS rank one update to a matrix. */
283 /* dcopy Level 1 BLAS that copies one vector to another . */
284 /* ddot Level 1 BLAS that computes the scalar product of two vectors. */
285 /* dnrm2 Level 1 BLAS that computes the norm of a vector. */
286 /* dscal Level 1 BLAS that scales a vector. */
287
288 /* \Remarks */
289
290 /* 1. Currently only HOWMNY = 'A' and 'P' are implemented. */
291
292 /* Let trans(X) denote the transpose of X. */
293
294 /* 2. Schur vectors are an orthogonal representation for the basis of */
295 /* Ritz vectors. Thus, their numerical properties are often superior. */
296 /* If RVEC = .TRUE. then the relationship */
297 /* A * V(:,1:IPARAM(5)) = V(:,1:IPARAM(5)) * T, and */
298 /* trans(V(:,1:IPARAM(5))) * V(:,1:IPARAM(5)) = I are approximately */
299 /* satisfied. Here T is the leading submatrix of order IPARAM(5) of the */
300 /* real upper quasi-triangular matrix stored workl(ipntr(12)). That is, */
301 /* T is block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; */
302 /* each 2-by-2 diagonal block has its diagonal elements equal and its */
303 /* off-diagonal elements of opposite sign. Corresponding to each 2-by-2 */
304 /* diagonal block is a complex conjugate pair of Ritz values. The real */
305 /* Ritz values are stored on the diagonal of T. */
306
307 /* 3. If IPARAM(7) = 3 or 4 and SIGMAI is not equal zero, then the user must */
308 /* form the IPARAM(5) Rayleigh quotients in order to transform the Ritz */
309 /* values computed by DNAUPD for OP to those of A*z = lambda*B*z. */
310 /* Set RVEC = .true. and HOWMNY = 'A', and */
311 /* compute */
312 /* trans(Z(:,I)) * A * Z(:,I) if DI(I) = 0. */
313 /* If DI(I) is not equal to zero and DI(I+1) = - D(I), */
314 /* then the desired real and imaginary parts of the Ritz value are */
315 /* trans(Z(:,I)) * A * Z(:,I) + trans(Z(:,I+1)) * A * Z(:,I+1), */
316 /* trans(Z(:,I)) * A * Z(:,I+1) - trans(Z(:,I+1)) * A * Z(:,I), */
317 /* respectively. */
318 /* Another possibility is to set RVEC = .true. and HOWMNY = 'P' and */
319 /* compute trans(V(:,1:IPARAM(5))) * A * V(:,1:IPARAM(5)) and then an upper */
320 /* quasi-triangular matrix of order IPARAM(5) is computed. See remark */
321 /* 2 above. */
322
323 /* \Authors */
324 /* Danny Sorensen Phuong Vu */
325 /* Richard Lehoucq CRPC / Rice University */
326 /* Chao Yang Houston, Texas */
327 /* Dept. of Computational & */
328 /* Applied Mathematics */
329 /* Rice University */
330 /* Houston, Texas */
331
332 /* \SCCS Information: @(#) */
333 /* FILE: neupd.F SID: 2.7 DATE OF SID: 09/20/00 RELEASE: 2 */
334
335 /* \EndLib */
336
337 /* ----------------------------------------------------------------------- */
dneupd_(logical * rvec,char * howmny,logical * select,doublereal * dr,doublereal * di,doublereal * z__,integer * ldz,doublereal * sigmar,doublereal * sigmai,doublereal * workev,char * bmat,integer * n,char * which,integer * nev,doublereal * tol,doublereal * resid,integer * ncv,doublereal * v,integer * ldv,integer * iparam,integer * ipntr,doublereal * workd,doublereal * workl,integer * lworkl,integer * info,ftnlen howmny_len,ftnlen bmat_len,ftnlen which_len)338 /* Subroutine */ int dneupd_(logical *rvec, char *howmny, logical *select,
339 doublereal *dr, doublereal *di, doublereal *z__, integer *ldz,
340 doublereal *sigmar, doublereal *sigmai, doublereal *workev, char *
341 bmat, integer *n, char *which, integer *nev, doublereal *tol,
342 doublereal *resid, integer *ncv, doublereal *v, integer *ldv, integer
343 *iparam, integer *ipntr, doublereal *workd, doublereal *workl,
344 integer *lworkl, integer *info, ftnlen howmny_len, ftnlen bmat_len,
345 ftnlen which_len)
346 {
347 /* System generated locals */
348 integer v_dim1, v_offset, z_dim1, z_offset, i__1;
349 doublereal d__1, d__2;
350
351 /* Builtin functions */
352 double pow_dd(doublereal *, doublereal *);
353 integer s_cmp(char *, char *, ftnlen, ftnlen);
354 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
355
356 /* Local variables */
357 static integer j, k, ih, jj, np;
358 static doublereal vl[1] /* was [1][1] */;
359 static integer ibd, ldh, ldq, iri;
360 static doublereal sep;
361 static integer irr, wri, wrr;
362 extern /* Subroutine */ int dger_(integer *, integer *, doublereal *,
363 doublereal *, integer *, doublereal *, integer *, doublereal *,
364 integer *);
365 static integer mode;
366 static doublereal eps23;
367 static integer ierr;
368 static doublereal temp;
369 static integer iwev;
370 static char type__[6];
371 extern doublereal dnrm2_(integer *, doublereal *, integer *);
372 static doublereal temp1;
373 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
374 integer *);
375 static integer ihbds, iconj;
376 extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
377 doublereal *, doublereal *, integer *, doublereal *, integer *,
378 doublereal *, doublereal *, integer *, ftnlen);
379 static doublereal conds;
380 static logical reord;
381 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
382 doublereal *, integer *);
383 static integer nconv;
384 extern /* Subroutine */ int dtrmm_(char *, char *, char *, char *,
385 integer *, integer *, doublereal *, doublereal *, integer *,
386 doublereal *, integer *, ftnlen, ftnlen, ftnlen, ftnlen), dmout_(
387 integer *, integer *, integer *, doublereal *, integer *, integer
388 *, char *, ftnlen);
389 static integer iwork[1];
390 static doublereal rnorm;
391 static integer ritzi;
392 extern /* Subroutine */ int dvout_(integer *, integer *, doublereal *,
393 integer *, char *, ftnlen), ivout_(integer *, integer *, integer *
394 , integer *, char *, ftnlen);
395 static integer ritzr;
396 extern /* Subroutine */ int dgeqr2_(integer *, integer *, doublereal *,
397 integer *, doublereal *, doublereal *, integer *);
398 extern doublereal dlapy2_(doublereal *, doublereal *);
399 extern /* Subroutine */ int dorm2r_(char *, char *, integer *, integer *,
400 integer *, doublereal *, integer *, doublereal *, doublereal *,
401 integer *, doublereal *, integer *, ftnlen, ftnlen);
402 extern doublereal dlamch_(char *, ftnlen);
403 static integer iheigi, iheigr, bounds, invsub, iuptri, msglvl, outncv,
404 ishift, numcnv;
405 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
406 doublereal *, integer *, doublereal *, integer *, ftnlen),
407 dlahqr_(logical *, logical *, integer *, integer *, integer *,
408 doublereal *, integer *, doublereal *, doublereal *, integer *,
409 integer *, doublereal *, integer *, integer *), dlaset_(char *,
410 integer *, integer *, doublereal *, doublereal *, doublereal *,
411 integer *, ftnlen), dtrevc_(char *, char *, logical *, integer *,
412 doublereal *, integer *, doublereal *, integer *, doublereal *,
413 integer *, integer *, integer *, doublereal *, integer *, ftnlen,
414 ftnlen), dtrsen_(char *, char *, logical *, integer *, doublereal
415 *, integer *, doublereal *, integer *, doublereal *, doublereal *,
416 integer *, doublereal *, doublereal *, doublereal *, integer *,
417 integer *, integer *, integer *, ftnlen, ftnlen), dngets_(integer
418 *, char *, integer *, integer *, doublereal *, doublereal *,
419 doublereal *, doublereal *, doublereal *, ftnlen);
420
421
422 /* %----------------------------------------------------% */
423 /* | Include files for debugging and timing information | */
424 /* %----------------------------------------------------% */
425
426
427 /* \SCCS Information: @(#) */
428 /* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
429
430 /* %---------------------------------% */
431 /* | See debug.doc for documentation | */
432 /* %---------------------------------% */
433
434 /* %------------------% */
435 /* | Scalar Arguments | */
436 /* %------------------% */
437
438 /* %--------------------------------% */
439 /* | See stat.doc for documentation | */
440 /* %--------------------------------% */
441
442 /* \SCCS Information: @(#) */
443 /* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
444
445
446
447 /* %-----------------% */
448 /* | Array Arguments | */
449 /* %-----------------% */
450
451
452 /* %------------% */
453 /* | Parameters | */
454 /* %------------% */
455
456
457 /* %---------------% */
458 /* | Local Scalars | */
459 /* %---------------% */
460
461
462 /* %----------------------% */
463 /* | External Subroutines | */
464 /* %----------------------% */
465
466
467 /* %--------------------% */
468 /* | External Functions | */
469 /* %--------------------% */
470
471
472 /* %---------------------% */
473 /* | Intrinsic Functions | */
474 /* %---------------------% */
475
476
477 /* %-----------------------% */
478 /* | Executable Statements | */
479 /* %-----------------------% */
480
481 /* %------------------------% */
482 /* | Set default parameters | */
483 /* %------------------------% */
484
485 /* Parameter adjustments */
486 z_dim1 = *ldz;
487 z_offset = 1 + z_dim1;
488 z__ -= z_offset;
489 --workd;
490 --resid;
491 --di;
492 --dr;
493 --workev;
494 --select;
495 v_dim1 = *ldv;
496 v_offset = 1 + v_dim1;
497 v -= v_offset;
498 --iparam;
499 --ipntr;
500 --workl;
501
502 /* Function Body */
503 msglvl = debug_1.mneupd;
504 mode = iparam[7];
505 nconv = iparam[5];
506 *info = 0;
507
508 /* %---------------------------------% */
509 /* | Get machine dependent constant. | */
510 /* %---------------------------------% */
511
512 eps23 = dlamch_("Epsilon-Machine", (ftnlen)15);
513 eps23 = pow_dd(&eps23, &c_b3);
514
515 /* %--------------% */
516 /* | Quick return | */
517 /* %--------------% */
518
519 ierr = 0;
520
521 if (nconv <= 0) {
522 ierr = -14;
523 } else if (*n <= 0) {
524 ierr = -1;
525 } else if (*nev <= 0) {
526 ierr = -2;
527 } else if (*ncv <= *nev + 1 || *ncv > *n) {
528 ierr = -3;
529 } else if (s_cmp(which, "LM", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which,
530 "SM", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which, "LR", (ftnlen)2,
531 (ftnlen)2) != 0 && s_cmp(which, "SR", (ftnlen)2, (ftnlen)2) != 0
532 && s_cmp(which, "LI", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which,
533 "SI", (ftnlen)2, (ftnlen)2) != 0) {
534 ierr = -5;
535 } else if (*(unsigned char *)bmat != 'I' && *(unsigned char *)bmat != 'G')
536 {
537 ierr = -6;
538 } else /* if(complicated condition) */ {
539 /* Computing 2nd power */
540 i__1 = *ncv;
541 if (*lworkl < i__1 * i__1 * 3 + *ncv * 6) {
542 ierr = -7;
543 } else if (*(unsigned char *)howmny != 'A' && *(unsigned char *)
544 howmny != 'P' && *(unsigned char *)howmny != 'S' && *rvec) {
545 ierr = -13;
546 } else if (*(unsigned char *)howmny == 'S') {
547 ierr = -12;
548 }
549 }
550
551 if (mode == 1 || mode == 2) {
552 s_copy(type__, "REGULR", (ftnlen)6, (ftnlen)6);
553 } else if (mode == 3 && *sigmai == 0.) {
554 s_copy(type__, "SHIFTI", (ftnlen)6, (ftnlen)6);
555 } else if (mode == 3) {
556 s_copy(type__, "REALPT", (ftnlen)6, (ftnlen)6);
557 } else if (mode == 4) {
558 s_copy(type__, "IMAGPT", (ftnlen)6, (ftnlen)6);
559 } else {
560 ierr = -10;
561 }
562 if (mode == 1 && *(unsigned char *)bmat == 'G') {
563 ierr = -11;
564 }
565
566 /* %------------% */
567 /* | Error Exit | */
568 /* %------------% */
569
570 if (ierr != 0) {
571 *info = ierr;
572 goto L9000;
573 }
574
575 /* %--------------------------------------------------------% */
576 /* | Pointer into WORKL for address of H, RITZ, BOUNDS, Q | */
577 /* | etc... and the remaining workspace. | */
578 /* | Also update pointer to be used on output. | */
579 /* | Memory is laid out as follows: | */
580 /* | workl(1:ncv*ncv) := generated Hessenberg matrix | */
581 /* | workl(ncv*ncv+1:ncv*ncv+2*ncv) := real and imaginary | */
582 /* | parts of ritz values | */
583 /* | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := error bounds | */
584 /* %--------------------------------------------------------% */
585
586 /* %-----------------------------------------------------------% */
587 /* | The following is used and set by DNEUPD . | */
588 /* | workl(ncv*ncv+3*ncv+1:ncv*ncv+4*ncv) := The untransformed | */
589 /* | real part of the Ritz values. | */
590 /* | workl(ncv*ncv+4*ncv+1:ncv*ncv+5*ncv) := The untransformed | */
591 /* | imaginary part of the Ritz values. | */
592 /* | workl(ncv*ncv+5*ncv+1:ncv*ncv+6*ncv) := The untransformed | */
593 /* | error bounds of the Ritz values | */
594 /* | workl(ncv*ncv+6*ncv+1:2*ncv*ncv+6*ncv) := Holds the upper | */
595 /* | quasi-triangular matrix for H | */
596 /* | workl(2*ncv*ncv+6*ncv+1: 3*ncv*ncv+6*ncv) := Holds the | */
597 /* | associated matrix representation of the invariant | */
598 /* | subspace for H. | */
599 /* | GRAND total of NCV * ( 3 * NCV + 6 ) locations. | */
600 /* %-----------------------------------------------------------% */
601
602 ih = ipntr[5];
603 ritzr = ipntr[6];
604 ritzi = ipntr[7];
605 bounds = ipntr[8];
606 ldh = *ncv;
607 ldq = *ncv;
608 iheigr = bounds + ldh;
609 iheigi = iheigr + ldh;
610 ihbds = iheigi + ldh;
611 iuptri = ihbds + ldh;
612 invsub = iuptri + ldh * *ncv;
613 ipntr[9] = iheigr;
614 ipntr[10] = iheigi;
615 ipntr[11] = ihbds;
616 ipntr[12] = iuptri;
617 ipntr[13] = invsub;
618 wrr = 1;
619 wri = *ncv + 1;
620 iwev = wri + *ncv;
621
622 /* %-----------------------------------------% */
623 /* | irr points to the REAL part of the Ritz | */
624 /* | values computed by _neigh before | */
625 /* | exiting _naup2. | */
626 /* | iri points to the IMAGINARY part of the | */
627 /* | Ritz values computed by _neigh | */
628 /* | before exiting _naup2. | */
629 /* | ibd points to the Ritz estimates | */
630 /* | computed by _neigh before exiting | */
631 /* | _naup2. | */
632 /* %-----------------------------------------% */
633
634 irr = ipntr[14] + *ncv * *ncv;
635 iri = irr + *ncv;
636 ibd = iri + *ncv;
637
638 /* %------------------------------------% */
639 /* | RNORM is B-norm of the RESID(1:N). | */
640 /* %------------------------------------% */
641
642 rnorm = workl[ih + 2];
643 workl[ih + 2] = 0.;
644
645 if (msglvl > 2) {
646 dvout_(&debug_1.logfil, ncv, &workl[irr], &debug_1.ndigit, "_neupd: "
647 "Real part of Ritz values passed in from _NAUPD.", (ftnlen)55);
648 dvout_(&debug_1.logfil, ncv, &workl[iri], &debug_1.ndigit, "_neupd: "
649 "Imag part of Ritz values passed in from _NAUPD.", (ftnlen)55);
650 dvout_(&debug_1.logfil, ncv, &workl[ibd], &debug_1.ndigit, "_neupd: "
651 "Ritz estimates passed in from _NAUPD.", (ftnlen)45);
652 }
653
654 if (*rvec) {
655
656 reord = FALSE_;
657
658 /* %---------------------------------------------------% */
659 /* | Use the temporary bounds array to store indices | */
660 /* | These will be used to mark the select array later | */
661 /* %---------------------------------------------------% */
662
663 i__1 = *ncv;
664 for (j = 1; j <= i__1; ++j) {
665 workl[bounds + j - 1] = (doublereal) j;
666 select[j] = FALSE_;
667 /* L10: */
668 }
669
670 /* %-------------------------------------% */
671 /* | Select the wanted Ritz values. | */
672 /* | Sort the Ritz values so that the | */
673 /* | wanted ones appear at the tailing | */
674 /* | NEV positions of workl(irr) and | */
675 /* | workl(iri). Move the corresponding | */
676 /* | error estimates in workl(bound) | */
677 /* | accordingly. | */
678 /* %-------------------------------------% */
679
680 np = *ncv - *nev;
681 ishift = 0;
682 dngets_(&ishift, which, nev, &np, &workl[irr], &workl[iri], &workl[
683 bounds], &workl[1], &workl[np + 1], (ftnlen)2);
684
685 if (msglvl > 2) {
686 dvout_(&debug_1.logfil, ncv, &workl[irr], &debug_1.ndigit, "_neu"
687 "pd: Real part of Ritz values after calling _NGETS.", (
688 ftnlen)54);
689 dvout_(&debug_1.logfil, ncv, &workl[iri], &debug_1.ndigit, "_neu"
690 "pd: Imag part of Ritz values after calling _NGETS.", (
691 ftnlen)54);
692 dvout_(&debug_1.logfil, ncv, &workl[bounds], &debug_1.ndigit,
693 "_neupd: Ritz value indices after calling _NGETS.", (
694 ftnlen)48);
695 }
696
697 /* %-----------------------------------------------------% */
698 /* | Record indices of the converged wanted Ritz values | */
699 /* | Mark the select array for possible reordering | */
700 /* %-----------------------------------------------------% */
701
702 numcnv = 0;
703 i__1 = *ncv;
704 for (j = 1; j <= i__1; ++j) {
705 /* Computing MAX */
706 d__1 = eps23, d__2 = dlapy2_(&workl[irr + *ncv - j], &workl[iri +
707 *ncv - j]);
708 temp1 = max(d__1,d__2);
709 jj = (integer) workl[bounds + *ncv - j];
710 if (numcnv < nconv && workl[ibd + jj - 1] <= *tol * temp1) {
711 select[jj] = TRUE_;
712 ++numcnv;
713 if (jj > *nev) {
714 reord = TRUE_;
715 }
716 }
717 /* L11: */
718 }
719
720 /* %-----------------------------------------------------------% */
721 /* | Check the count (numcnv) of converged Ritz values with | */
722 /* | the number (nconv) reported by dnaupd. If these two | */
723 /* | are different then there has probably been an error | */
724 /* | caused by incorrect passing of the dnaupd data. | */
725 /* %-----------------------------------------------------------% */
726
727 if (msglvl > 2) {
728 ivout_(&debug_1.logfil, &c__1, &numcnv, &debug_1.ndigit, "_neupd"
729 ": Number of specified eigenvalues", (ftnlen)39);
730 ivout_(&debug_1.logfil, &c__1, &nconv, &debug_1.ndigit, "_neupd:"
731 " Number of \"converged\" eigenvalues", (ftnlen)41);
732 }
733
734 if (numcnv != nconv) {
735 *info = -15;
736 goto L9000;
737 }
738
739 /* %-----------------------------------------------------------% */
740 /* | Call LAPACK routine dlahqr to compute the real Schur form | */
741 /* | of the upper Hessenberg matrix returned by DNAUPD . | */
742 /* | Make a copy of the upper Hessenberg matrix. | */
743 /* | Initialize the Schur vector matrix Q to the identity. | */
744 /* %-----------------------------------------------------------% */
745
746 i__1 = ldh * *ncv;
747 dcopy_(&i__1, &workl[ih], &c__1, &workl[iuptri], &c__1);
748 dlaset_("All", ncv, ncv, &c_b37, &c_b38, &workl[invsub], &ldq, (
749 ftnlen)3);
750 dlahqr_(&c_true, &c_true, ncv, &c__1, ncv, &workl[iuptri], &ldh, &
751 workl[iheigr], &workl[iheigi], &c__1, ncv, &workl[invsub], &
752 ldq, &ierr);
753 dcopy_(ncv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds], &c__1);
754
755 if (ierr != 0) {
756 *info = -8;
757 goto L9000;
758 }
759
760 if (msglvl > 1) {
761 dvout_(&debug_1.logfil, ncv, &workl[iheigr], &debug_1.ndigit,
762 "_neupd: Real part of the eigenvalues of H", (ftnlen)41);
763 dvout_(&debug_1.logfil, ncv, &workl[iheigi], &debug_1.ndigit,
764 "_neupd: Imaginary part of the Eigenvalues of H", (ftnlen)
765 46);
766 dvout_(&debug_1.logfil, ncv, &workl[ihbds], &debug_1.ndigit,
767 "_neupd: Last row of the Schur vector matrix", (ftnlen)43)
768 ;
769 if (msglvl > 3) {
770 dmout_(&debug_1.logfil, ncv, ncv, &workl[iuptri], &ldh, &
771 debug_1.ndigit, "_neupd: The upper quasi-triangular "
772 "matrix ", (ftnlen)42);
773 }
774 }
775
776 if (reord) {
777
778 /* %-----------------------------------------------------% */
779 /* | Reorder the computed upper quasi-triangular matrix. | */
780 /* %-----------------------------------------------------% */
781
782 dtrsen_("None", "V", &select[1], ncv, &workl[iuptri], &ldh, &
783 workl[invsub], &ldq, &workl[iheigr], &workl[iheigi], &
784 nconv, &conds, &sep, &workl[ihbds], ncv, iwork, &c__1, &
785 ierr, (ftnlen)4, (ftnlen)1);
786
787 if (ierr == 1) {
788 *info = 1;
789 goto L9000;
790 }
791
792 if (msglvl > 2) {
793 dvout_(&debug_1.logfil, ncv, &workl[iheigr], &debug_1.ndigit,
794 "_neupd: Real part of the eigenvalues of H--reordered"
795 , (ftnlen)52);
796 dvout_(&debug_1.logfil, ncv, &workl[iheigi], &debug_1.ndigit,
797 "_neupd: Imag part of the eigenvalues of H--reordered"
798 , (ftnlen)52);
799 if (msglvl > 3) {
800 dmout_(&debug_1.logfil, ncv, ncv, &workl[iuptri], &ldq, &
801 debug_1.ndigit, "_neupd: Quasi-triangular matrix"
802 " after re-ordering", (ftnlen)49);
803 }
804 }
805
806 }
807
808 /* %---------------------------------------% */
809 /* | Copy the last row of the Schur vector | */
810 /* | into workl(ihbds). This will be used | */
811 /* | to compute the Ritz estimates of | */
812 /* | converged Ritz values. | */
813 /* %---------------------------------------% */
814
815 dcopy_(ncv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds], &c__1);
816
817 /* %----------------------------------------------------% */
818 /* | Place the computed eigenvalues of H into DR and DI | */
819 /* | if a spectral transformation was not used. | */
820 /* %----------------------------------------------------% */
821
822 if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) == 0) {
823 dcopy_(&nconv, &workl[iheigr], &c__1, &dr[1], &c__1);
824 dcopy_(&nconv, &workl[iheigi], &c__1, &di[1], &c__1);
825 }
826
827 /* %----------------------------------------------------------% */
828 /* | Compute the QR factorization of the matrix representing | */
829 /* | the wanted invariant subspace located in the first NCONV | */
830 /* | columns of workl(invsub,ldq). | */
831 /* %----------------------------------------------------------% */
832
833 dgeqr2_(ncv, &nconv, &workl[invsub], &ldq, &workev[1], &workev[*ncv +
834 1], &ierr);
835
836 /* %---------------------------------------------------------% */
837 /* | * Postmultiply V by Q using dorm2r . | */
838 /* | * Copy the first NCONV columns of VQ into Z. | */
839 /* | * Postmultiply Z by R. | */
840 /* | The N by NCONV matrix Z is now a matrix representation | */
841 /* | of the approximate invariant subspace associated with | */
842 /* | the Ritz values in workl(iheigr) and workl(iheigi) | */
843 /* | The first NCONV columns of V are now approximate Schur | */
844 /* | vectors associated with the real upper quasi-triangular | */
845 /* | matrix of order NCONV in workl(iuptri) | */
846 /* %---------------------------------------------------------% */
847
848 dorm2r_("Right", "Notranspose", n, ncv, &nconv, &workl[invsub], &ldq,
849 &workev[1], &v[v_offset], ldv, &workd[*n + 1], &ierr, (ftnlen)
850 5, (ftnlen)11);
851 dlacpy_("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz, (
852 ftnlen)3);
853
854 i__1 = nconv;
855 for (j = 1; j <= i__1; ++j) {
856
857 /* %---------------------------------------------------% */
858 /* | Perform both a column and row scaling if the | */
859 /* | diagonal element of workl(invsub,ldq) is negative | */
860 /* | I'm lazy and don't take advantage of the upper | */
861 /* | quasi-triangular form of workl(iuptri,ldq) | */
862 /* | Note that since Q is orthogonal, R is a diagonal | */
863 /* | matrix consisting of plus or minus ones | */
864 /* %---------------------------------------------------% */
865
866 if (workl[invsub + (j - 1) * ldq + j - 1] < 0.) {
867 dscal_(&nconv, &c_b64, &workl[iuptri + j - 1], &ldq);
868 dscal_(&nconv, &c_b64, &workl[iuptri + (j - 1) * ldq], &c__1);
869 }
870
871 /* L20: */
872 }
873
874 if (*(unsigned char *)howmny == 'A') {
875
876 /* %--------------------------------------------% */
877 /* | Compute the NCONV wanted eigenvectors of T | */
878 /* | located in workl(iuptri,ldq). | */
879 /* %--------------------------------------------% */
880
881 i__1 = *ncv;
882 for (j = 1; j <= i__1; ++j) {
883 if (j <= nconv) {
884 select[j] = TRUE_;
885 } else {
886 select[j] = FALSE_;
887 }
888 /* L30: */
889 }
890
891 dtrevc_("Right", "Select", &select[1], ncv, &workl[iuptri], &ldq,
892 vl, &c__1, &workl[invsub], &ldq, ncv, &outncv, &workev[1],
893 &ierr, (ftnlen)5, (ftnlen)6);
894
895 if (ierr != 0) {
896 *info = -9;
897 goto L9000;
898 }
899
900 /* %------------------------------------------------% */
901 /* | Scale the returning eigenvectors so that their | */
902 /* | Euclidean norms are all one. LAPACK subroutine | */
903 /* | dtrevc returns each eigenvector normalized so | */
904 /* | that the element of largest magnitude has | */
905 /* | magnitude 1; | */
906 /* %------------------------------------------------% */
907
908 iconj = 0;
909 i__1 = nconv;
910 for (j = 1; j <= i__1; ++j) {
911
912 if (workl[iheigi + j - 1] == 0.) {
913
914 /* %----------------------% */
915 /* | real eigenvalue case | */
916 /* %----------------------% */
917
918 temp = dnrm2_(ncv, &workl[invsub + (j - 1) * ldq], &c__1);
919 d__1 = 1. / temp;
920 dscal_(ncv, &d__1, &workl[invsub + (j - 1) * ldq], &c__1);
921
922 } else {
923
924 /* %-------------------------------------------% */
925 /* | Complex conjugate pair case. Note that | */
926 /* | since the real and imaginary part of | */
927 /* | the eigenvector are stored in consecutive | */
928 /* | columns, we further normalize by the | */
929 /* | square root of two. | */
930 /* %-------------------------------------------% */
931
932 if (iconj == 0) {
933 d__1 = dnrm2_(ncv, &workl[invsub + (j - 1) * ldq], &
934 c__1);
935 d__2 = dnrm2_(ncv, &workl[invsub + j * ldq], &c__1);
936 temp = dlapy2_(&d__1, &d__2);
937 d__1 = 1. / temp;
938 dscal_(ncv, &d__1, &workl[invsub + (j - 1) * ldq], &
939 c__1);
940 d__1 = 1. / temp;
941 dscal_(ncv, &d__1, &workl[invsub + j * ldq], &c__1);
942 iconj = 1;
943 } else {
944 iconj = 0;
945 }
946
947 }
948
949 /* L40: */
950 }
951
952 dgemv_("T", ncv, &nconv, &c_b38, &workl[invsub], &ldq, &workl[
953 ihbds], &c__1, &c_b37, &workev[1], &c__1, (ftnlen)1);
954
955 iconj = 0;
956 i__1 = nconv;
957 for (j = 1; j <= i__1; ++j) {
958 if (workl[iheigi + j - 1] != 0.) {
959
960 /* %-------------------------------------------% */
961 /* | Complex conjugate pair case. Note that | */
962 /* | since the real and imaginary part of | */
963 /* | the eigenvector are stored in consecutive | */
964 /* %-------------------------------------------% */
965
966 if (iconj == 0) {
967 workev[j] = dlapy2_(&workev[j], &workev[j + 1]);
968 workev[j + 1] = workev[j];
969 iconj = 1;
970 } else {
971 iconj = 0;
972 }
973 }
974 /* L45: */
975 }
976
977 if (msglvl > 2) {
978 dcopy_(ncv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds], &
979 c__1);
980 dvout_(&debug_1.logfil, ncv, &workl[ihbds], &debug_1.ndigit,
981 "_neupd: Last row of the eigenvector matrix for T", (
982 ftnlen)48);
983 if (msglvl > 3) {
984 dmout_(&debug_1.logfil, ncv, ncv, &workl[invsub], &ldq, &
985 debug_1.ndigit, "_neupd: The eigenvector matrix "
986 "for T", (ftnlen)36);
987 }
988 }
989
990 /* %---------------------------------------% */
991 /* | Copy Ritz estimates into workl(ihbds) | */
992 /* %---------------------------------------% */
993
994 dcopy_(&nconv, &workev[1], &c__1, &workl[ihbds], &c__1);
995
996 /* %---------------------------------------------------------% */
997 /* | Compute the QR factorization of the eigenvector matrix | */
998 /* | associated with leading portion of T in the first NCONV | */
999 /* | columns of workl(invsub,ldq). | */
1000 /* %---------------------------------------------------------% */
1001
1002 dgeqr2_(ncv, &nconv, &workl[invsub], &ldq, &workev[1], &workev[*
1003 ncv + 1], &ierr);
1004
1005 /* %----------------------------------------------% */
1006 /* | * Postmultiply Z by Q. | */
1007 /* | * Postmultiply Z by R. | */
1008 /* | The N by NCONV matrix Z is now contains the | */
1009 /* | Ritz vectors associated with the Ritz values | */
1010 /* | in workl(iheigr) and workl(iheigi). | */
1011 /* %----------------------------------------------% */
1012
1013 dorm2r_("Right", "Notranspose", n, ncv, &nconv, &workl[invsub], &
1014 ldq, &workev[1], &z__[z_offset], ldz, &workd[*n + 1], &
1015 ierr, (ftnlen)5, (ftnlen)11);
1016
1017 dtrmm_("Right", "Upper", "No transpose", "Non-unit", n, &nconv, &
1018 c_b38, &workl[invsub], &ldq, &z__[z_offset], ldz, (ftnlen)
1019 5, (ftnlen)5, (ftnlen)12, (ftnlen)8);
1020
1021 }
1022
1023 } else {
1024
1025 /* %------------------------------------------------------% */
1026 /* | An approximate invariant subspace is not needed. | */
1027 /* | Place the Ritz values computed DNAUPD into DR and DI | */
1028 /* %------------------------------------------------------% */
1029
1030 dcopy_(&nconv, &workl[ritzr], &c__1, &dr[1], &c__1);
1031 dcopy_(&nconv, &workl[ritzi], &c__1, &di[1], &c__1);
1032 dcopy_(&nconv, &workl[ritzr], &c__1, &workl[iheigr], &c__1);
1033 dcopy_(&nconv, &workl[ritzi], &c__1, &workl[iheigi], &c__1);
1034 dcopy_(&nconv, &workl[bounds], &c__1, &workl[ihbds], &c__1);
1035 }
1036
1037 /* %------------------------------------------------% */
1038 /* | Transform the Ritz values and possibly vectors | */
1039 /* | and corresponding error bounds of OP to those | */
1040 /* | of A*x = lambda*B*x. | */
1041 /* %------------------------------------------------% */
1042
1043 if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) == 0) {
1044
1045 if (*rvec) {
1046 dscal_(ncv, &rnorm, &workl[ihbds], &c__1);
1047 }
1048
1049 } else {
1050
1051 /* %---------------------------------------% */
1052 /* | A spectral transformation was used. | */
1053 /* | * Determine the Ritz estimates of the | */
1054 /* | Ritz values in the original system. | */
1055 /* %---------------------------------------% */
1056
1057 if (s_cmp(type__, "SHIFTI", (ftnlen)6, (ftnlen)6) == 0) {
1058
1059 if (*rvec) {
1060 dscal_(ncv, &rnorm, &workl[ihbds], &c__1);
1061 }
1062
1063 i__1 = *ncv;
1064 for (k = 1; k <= i__1; ++k) {
1065 temp = dlapy2_(&workl[iheigr + k - 1], &workl[iheigi + k - 1])
1066 ;
1067 workl[ihbds + k - 1] = (d__1 = workl[ihbds + k - 1], abs(d__1)
1068 ) / temp / temp;
1069 /* L50: */
1070 }
1071
1072 } else if (s_cmp(type__, "REALPT", (ftnlen)6, (ftnlen)6) == 0) {
1073
1074 i__1 = *ncv;
1075 for (k = 1; k <= i__1; ++k) {
1076 /* L60: */
1077 }
1078
1079 } else if (s_cmp(type__, "IMAGPT", (ftnlen)6, (ftnlen)6) == 0) {
1080
1081 i__1 = *ncv;
1082 for (k = 1; k <= i__1; ++k) {
1083 /* L70: */
1084 }
1085
1086 }
1087
1088 /* %-----------------------------------------------------------% */
1089 /* | * Transform the Ritz values back to the original system. | */
1090 /* | For TYPE = 'SHIFTI' the transformation is | */
1091 /* | lambda = 1/theta + sigma | */
1092 /* | For TYPE = 'REALPT' or 'IMAGPT' the user must from | */
1093 /* | Rayleigh quotients or a projection. See remark 3 above.| */
1094 /* | NOTES: | */
1095 /* | *The Ritz vectors are not affected by the transformation. | */
1096 /* %-----------------------------------------------------------% */
1097
1098 if (s_cmp(type__, "SHIFTI", (ftnlen)6, (ftnlen)6) == 0) {
1099
1100 i__1 = *ncv;
1101 for (k = 1; k <= i__1; ++k) {
1102 temp = dlapy2_(&workl[iheigr + k - 1], &workl[iheigi + k - 1])
1103 ;
1104 workl[iheigr + k - 1] = workl[iheigr + k - 1] / temp / temp +
1105 *sigmar;
1106 workl[iheigi + k - 1] = -workl[iheigi + k - 1] / temp / temp
1107 + *sigmai;
1108 /* L80: */
1109 }
1110
1111 dcopy_(&nconv, &workl[iheigr], &c__1, &dr[1], &c__1);
1112 dcopy_(&nconv, &workl[iheigi], &c__1, &di[1], &c__1);
1113
1114 } else if (s_cmp(type__, "REALPT", (ftnlen)6, (ftnlen)6) == 0 ||
1115 s_cmp(type__, "IMAGPT", (ftnlen)6, (ftnlen)6) == 0) {
1116
1117 dcopy_(&nconv, &workl[iheigr], &c__1, &dr[1], &c__1);
1118 dcopy_(&nconv, &workl[iheigi], &c__1, &di[1], &c__1);
1119
1120 }
1121
1122 }
1123
1124 if (s_cmp(type__, "SHIFTI", (ftnlen)6, (ftnlen)6) == 0 && msglvl > 1) {
1125 dvout_(&debug_1.logfil, &nconv, &dr[1], &debug_1.ndigit, "_neupd: Un"
1126 "transformed real part of the Ritz valuess.", (ftnlen)52);
1127 dvout_(&debug_1.logfil, &nconv, &di[1], &debug_1.ndigit, "_neupd: Un"
1128 "transformed imag part of the Ritz valuess.", (ftnlen)52);
1129 dvout_(&debug_1.logfil, &nconv, &workl[ihbds], &debug_1.ndigit, "_ne"
1130 "upd: Ritz estimates of untransformed Ritz values.", (ftnlen)
1131 52);
1132 } else if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) == 0 && msglvl >
1133 1) {
1134 dvout_(&debug_1.logfil, &nconv, &dr[1], &debug_1.ndigit, "_neupd: Re"
1135 "al parts of converged Ritz values.", (ftnlen)44);
1136 dvout_(&debug_1.logfil, &nconv, &di[1], &debug_1.ndigit, "_neupd: Im"
1137 "ag parts of converged Ritz values.", (ftnlen)44);
1138 dvout_(&debug_1.logfil, &nconv, &workl[ihbds], &debug_1.ndigit, "_ne"
1139 "upd: Associated Ritz estimates.", (ftnlen)34);
1140 }
1141
1142 /* %-------------------------------------------------% */
1143 /* | Eigenvector Purification step. Formally perform | */
1144 /* | one of inverse subspace iteration. Only used | */
1145 /* | for MODE = 2. | */
1146 /* %-------------------------------------------------% */
1147
1148 if (*rvec && *(unsigned char *)howmny == 'A' && s_cmp(type__, "SHIFTI", (
1149 ftnlen)6, (ftnlen)6) == 0) {
1150
1151 /* %------------------------------------------------% */
1152 /* | Purify the computed Ritz vectors by adding a | */
1153 /* | little bit of the residual vector: | */
1154 /* | T | */
1155 /* | resid(:)*( e s ) / theta | */
1156 /* | NCV | */
1157 /* | where H s = s theta. Remember that when theta | */
1158 /* | has nonzero imaginary part, the corresponding | */
1159 /* | Ritz vector is stored across two columns of Z. | */
1160 /* %------------------------------------------------% */
1161
1162 iconj = 0;
1163 i__1 = nconv;
1164 for (j = 1; j <= i__1; ++j) {
1165 if (workl[iheigi + j - 1] == 0.) {
1166 workev[j] = workl[invsub + (j - 1) * ldq + *ncv - 1] / workl[
1167 iheigr + j - 1];
1168 } else if (iconj == 0) {
1169 temp = dlapy2_(&workl[iheigr + j - 1], &workl[iheigi + j - 1])
1170 ;
1171 workev[j] = (workl[invsub + (j - 1) * ldq + *ncv - 1] * workl[
1172 iheigr + j - 1] + workl[invsub + j * ldq + *ncv - 1] *
1173 workl[iheigi + j - 1]) / temp / temp;
1174 workev[j + 1] = (workl[invsub + j * ldq + *ncv - 1] * workl[
1175 iheigr + j - 1] - workl[invsub + (j - 1) * ldq + *ncv
1176 - 1] * workl[iheigi + j - 1]) / temp / temp;
1177 iconj = 1;
1178 } else {
1179 iconj = 0;
1180 }
1181 /* L110: */
1182 }
1183
1184 /* %---------------------------------------% */
1185 /* | Perform a rank one update to Z and | */
1186 /* | purify all the Ritz vectors together. | */
1187 /* %---------------------------------------% */
1188
1189 dger_(n, &nconv, &c_b38, &resid[1], &c__1, &workev[1], &c__1, &z__[
1190 z_offset], ldz);
1191
1192 }
1193
1194 L9000:
1195
1196 return 0;
1197
1198 /* %---------------% */
1199 /* | End of DNEUPD | */
1200 /* %---------------% */
1201
1202 } /* dneupd_ */
1203
1204