1 /* ./src_f77/zneupd.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_ZNEUPD;
15
16 #define debug_1 debug_ZNEUPD
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_ZNEUPD;
25
26 #define timing_1 timing_ZNEUPD
27
28 /* Table of constant values */
29
30 static doublecomplex c_b1 = {1.,0.};
31 static doublecomplex c_b2 = {0.,0.};
32 static doublereal c_b5 = .66666666666666663;
33 static integer c__1 = 1;
34 static logical c_true = TRUE_;
35
36 /* \BeginDoc */
37
38 /* \Name: zneupd */
39
40 /* \Description: */
41 /* This subroutine returns the converged approximations to eigenvalues */
42 /* of A*z = lambda*B*z and (optionally): */
43
44 /* (1) The corresponding approximate eigenvectors; */
45
46 /* (2) An orthonormal basis for the associated approximate */
47 /* invariant subspace; */
48
49 /* (3) Both. */
50
51 /* There is negligible additional cost to obtain eigenvectors. An orthonormal */
52 /* basis is always computed. There is an additional storage cost of n*nev */
53 /* if both are requested (in this case a separate array Z must be supplied). */
54
55 /* The approximate eigenvalues and eigenvectors of A*z = lambda*B*z */
56 /* are derived from approximate eigenvalues and eigenvectors of */
57 /* of the linear operator OP prescribed by the MODE selection in the */
58 /* call to ZNAUPD. ZNAUPD must be called before this routine is called. */
59 /* These approximate eigenvalues and vectors are commonly called Ritz */
60 /* values and Ritz vectors respectively. They are referred to as such */
61 /* in the comments that follow. The computed orthonormal basis for the */
62 /* invariant subspace corresponding to these Ritz values is referred to as a */
63 /* Schur basis. */
64
65 /* The definition of OP as well as other terms and the relation of computed */
66 /* Ritz values and vectors of OP with respect to the given problem */
67 /* A*z = lambda*B*z may be found in the header of ZNAUPD. For a brief */
68 /* description, see definitions of IPARAM(7), MODE and WHICH in the */
69 /* documentation of ZNAUPD. */
70
71 /* \Usage: */
72 /* call zneupd */
73 /* ( RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, WORKEV, BMAT, */
74 /* N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, */
75 /* WORKL, LWORKL, RWORK, INFO ) */
76
77 /* \Arguments: */
78 /* RVEC LOGICAL (INPUT) */
79 /* Specifies whether a basis for the invariant subspace corresponding */
80 /* to the converged Ritz value approximations for the eigenproblem */
81 /* A*z = lambda*B*z is computed. */
82
83 /* RVEC = .FALSE. Compute Ritz values only. */
84
85 /* RVEC = .TRUE. Compute Ritz vectors or Schur vectors. */
86 /* See Remarks below. */
87
88 /* HOWMNY Character*1 (INPUT) */
89 /* Specifies the form of the basis for the invariant subspace */
90 /* corresponding to the converged Ritz values that is to be computed. */
91
92 /* = 'A': Compute NEV Ritz vectors; */
93 /* = 'P': Compute NEV Schur vectors; */
94 /* = 'S': compute some of the Ritz vectors, specified */
95 /* by the logical array SELECT. */
96
97 /* SELECT Logical array of dimension NCV. (INPUT) */
98 /* If HOWMNY = 'S', SELECT specifies the Ritz vectors to be */
99 /* computed. To select the Ritz vector corresponding to a */
100 /* Ritz value D(j), SELECT(j) must be set to .TRUE.. */
101 /* If HOWMNY = 'A' or 'P', SELECT need not be initialized */
102 /* but it is used as internal workspace. */
103
104 /* D Complex*16 array of dimension NEV+1. (OUTPUT) */
105 /* On exit, D contains the Ritz approximations */
106 /* to the eigenvalues lambda for A*z = lambda*B*z. */
107
108 /* Z Complex*16 N by NEV array (OUTPUT) */
109 /* On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of */
110 /* Z represents approximate eigenvectors (Ritz vectors) corresponding */
111 /* to the NCONV=IPARAM(5) Ritz values for eigensystem */
112 /* A*z = lambda*B*z. */
113
114 /* If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED. */
115
116 /* NOTE: If if RVEC = .TRUE. and a Schur basis is not required, */
117 /* the array Z may be set equal to first NEV+1 columns of the Arnoldi */
118 /* basis array V computed by ZNAUPD. In this case the Arnoldi basis */
119 /* will be destroyed and overwritten with the eigenvector basis. */
120
121 /* LDZ Integer. (INPUT) */
122 /* The leading dimension of the array Z. If Ritz vectors are */
123 /* desired, then LDZ .ge. max( 1, N ) is required. */
124 /* In any case, LDZ .ge. 1 is required. */
125
126 /* SIGMA Complex*16 (INPUT) */
127 /* If IPARAM(7) = 3 then SIGMA represents the shift. */
128 /* Not referenced if IPARAM(7) = 1 or 2. */
129
130 /* WORKEV Complex*16 work array of dimension 2*NCV. (WORKSPACE) */
131
132 /* **** The remaining arguments MUST be the same as for the **** */
133 /* **** call to ZNAUPD that was just completed. **** */
134
135 /* NOTE: The remaining arguments */
136
137 /* BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, */
138 /* WORKD, WORKL, LWORKL, RWORK, INFO */
139
140 /* must be passed directly to ZNEUPD following the last call */
141 /* to ZNAUPD. These arguments MUST NOT BE MODIFIED between */
142 /* the the last call to ZNAUPD and the call to ZNEUPD. */
143
144 /* Three of these parameters (V, WORKL and INFO) are also output parameters: */
145
146 /* V Complex*16 N by NCV array. (INPUT/OUTPUT) */
147
148 /* Upon INPUT: the NCV columns of V contain the Arnoldi basis */
149 /* vectors for OP as constructed by ZNAUPD . */
150
151 /* Upon OUTPUT: If RVEC = .TRUE. the first NCONV=IPARAM(5) columns */
152 /* contain approximate Schur vectors that span the */
153 /* desired invariant subspace. */
154
155 /* NOTE: If the array Z has been set equal to first NEV+1 columns */
156 /* of the array V and RVEC=.TRUE. and HOWMNY= 'A', then the */
157 /* Arnoldi basis held by V has been overwritten by the desired */
158 /* Ritz vectors. If a separate array Z has been passed then */
159 /* the first NCONV=IPARAM(5) columns of V will contain approximate */
160 /* Schur vectors that span the desired invariant subspace. */
161
162 /* WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE) */
163 /* WORKL(1:ncv*ncv+2*ncv) contains information obtained in */
164 /* znaupd. They are not changed by zneupd. */
165 /* WORKL(ncv*ncv+2*ncv+1:3*ncv*ncv+4*ncv) holds the */
166 /* untransformed Ritz values, the untransformed error estimates of */
167 /* the Ritz values, the upper triangular matrix for H, and the */
168 /* associated matrix representation of the invariant subspace for H. */
169
170 /* Note: IPNTR(9:13) contains the pointer into WORKL for addresses */
171 /* of the above information computed by zneupd. */
172 /* ------------------------------------------------------------- */
173 /* IPNTR(9): pointer to the NCV RITZ values of the */
174 /* original system. */
175 /* IPNTR(10): Not used */
176 /* IPNTR(11): pointer to the NCV corresponding error estimates. */
177 /* IPNTR(12): pointer to the NCV by NCV upper triangular */
178 /* Schur matrix for H. */
179 /* IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors */
180 /* of the upper Hessenberg matrix H. Only referenced by */
181 /* zneupd if RVEC = .TRUE. See Remark 2 below. */
182 /* ------------------------------------------------------------- */
183
184 /* INFO Integer. (OUTPUT) */
185 /* Error flag on output. */
186 /* = 0: Normal exit. */
187
188 /* = 1: The Schur form computed by LAPACK routine csheqr */
189 /* could not be reordered by LAPACK routine ztrsen. */
190 /* Re-enter subroutine zneupd with IPARAM(5)=NCV and */
191 /* increase the size of the array D to have */
192 /* dimension at least dimension NCV and allocate at least NCV */
193 /* columns for Z. NOTE: Not necessary if Z and V share */
194 /* the same space. Please notify the authors if this error */
195 /* occurs. */
196
197 /* = -1: N must be positive. */
198 /* = -2: NEV must be positive. */
199 /* = -3: NCV-NEV >= 1 and less than or equal to N. */
200 /* = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI' */
201 /* = -6: BMAT must be one of 'I' or 'G'. */
202 /* = -7: Length of private work WORKL array is not sufficient. */
203 /* = -8: Error return from LAPACK eigenvalue calculation. */
204 /* This should never happened. */
205 /* = -9: Error return from calculation of eigenvectors. */
206 /* Informational error from LAPACK routine ztrevc. */
207 /* = -10: IPARAM(7) must be 1,2,3 */
208 /* = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible. */
209 /* = -12: HOWMNY = 'S' not yet implemented */
210 /* = -13: HOWMNY must be one of 'A' or 'P' if RVEC = .true. */
211 /* = -14: ZNAUPD did not find any eigenvalues to sufficient */
212 /* accuracy. */
213 /* = -15: ZNEUPD got a different count of the number of converged */
214 /* Ritz values than ZNAUPD got. This indicates the user */
215 /* probably made an error in passing data from ZNAUPD to */
216 /* ZNEUPD or that the data was modified before entering */
217 /* ZNEUPD */
218
219 /* \BeginLib */
220
221 /* \References: */
222 /* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
223 /* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
224 /* pp 357-385. */
225 /* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
226 /* Restarted Arnoldi Iteration", Rice University Technical Report */
227 /* TR95-13, Department of Computational and Applied Mathematics. */
228 /* 3. B. Nour-Omid, B. N. Parlett, T. Ericsson and P. S. Jensen, */
229 /* "How to Implement the Spectral Transformation", Math Comp., */
230 /* Vol. 48, No. 178, April, 1987 pp. 664-673. */
231
232 /* \Routines called: */
233 /* ivout ARPACK utility routine that prints integers. */
234 /* zmout ARPACK utility routine that prints matrices */
235 /* zvout ARPACK utility routine that prints vectors. */
236 /* zgeqr2 LAPACK routine that computes the QR factorization of */
237 /* a matrix. */
238 /* zlacpy LAPACK matrix copy routine. */
239 /* zlahqr LAPACK routine that computes the Schur form of a */
240 /* upper Hessenberg matrix. */
241 /* zlaset LAPACK matrix initialization routine. */
242 /* ztrevc LAPACK routine to compute the eigenvectors of a matrix */
243 /* in upper triangular form. */
244 /* ztrsen LAPACK routine that re-orders the Schur form. */
245 /* zunm2r LAPACK routine that applies an orthogonal matrix in */
246 /* factored form. */
247 /* dlamch LAPACK routine that determines machine constants. */
248 /* ztrmm Level 3 BLAS matrix times an upper triangular matrix. */
249 /* zgeru Level 2 BLAS rank one update to a matrix. */
250 /* zcopy Level 1 BLAS that copies one vector to another . */
251 /* zscal Level 1 BLAS that scales a vector. */
252 /* zdscal Level 1 BLAS that scales a complex vector by a real number. */
253 /* dznrm2 Level 1 BLAS that computes the norm of a complex vector. */
254
255 /* \Remarks */
256
257 /* 1. Currently only HOWMNY = 'A' and 'P' are implemented. */
258
259 /* 2. Schur vectors are an orthogonal representation for the basis of */
260 /* Ritz vectors. Thus, their numerical properties are often superior. */
261 /* If RVEC = .true. then the relationship */
262 /* A * V(:,1:IPARAM(5)) = V(:,1:IPARAM(5)) * T, and */
263 /* transpose( V(:,1:IPARAM(5)) ) * V(:,1:IPARAM(5)) = I */
264 /* are approximately satisfied. */
265 /* Here T is the leading submatrix of order IPARAM(5) of the */
266 /* upper triangular matrix stored workl(ipntr(12)). */
267
268 /* \Authors */
269 /* Danny Sorensen Phuong Vu */
270 /* Richard Lehoucq CRPC / Rice University */
271 /* Chao Yang Houston, Texas */
272 /* Dept. of Computational & */
273 /* Applied Mathematics */
274 /* Rice University */
275 /* Houston, Texas */
276
277 /* \SCCS Information: @(#) */
278 /* FILE: neupd.F SID: 2.8 DATE OF SID: 07/21/02 RELEASE: 2 */
279
280 /* \EndLib */
281
282 /* ----------------------------------------------------------------------- */
zneupd_(logical * rvec,char * howmny,logical * select,doublecomplex * d__,doublecomplex * z__,integer * ldz,doublecomplex * sigma,doublecomplex * workev,char * bmat,integer * n,char * which,integer * nev,doublereal * tol,doublecomplex * resid,integer * ncv,doublecomplex * v,integer * ldv,integer * iparam,integer * ipntr,doublecomplex * workd,doublecomplex * workl,integer * lworkl,doublereal * rwork,integer * info,ftnlen howmny_len,ftnlen bmat_len,ftnlen which_len)283 /* Subroutine */ int zneupd_(logical *rvec, char *howmny, logical *select,
284 doublecomplex *d__, doublecomplex *z__, integer *ldz, doublecomplex *
285 sigma, doublecomplex *workev, char *bmat, integer *n, char *which,
286 integer *nev, doublereal *tol, doublecomplex *resid, integer *ncv,
287 doublecomplex *v, integer *ldv, integer *iparam, integer *ipntr,
288 doublecomplex *workd, doublecomplex *workl, integer *lworkl,
289 doublereal *rwork, integer *info, ftnlen howmny_len, ftnlen bmat_len,
290 ftnlen which_len)
291 {
292 /* System generated locals */
293 integer v_dim1, v_offset, z_dim1, z_offset, i__1, i__2;
294 doublereal d__1, d__2, d__3, d__4;
295 doublecomplex z__1, z__2;
296
297 /* Builtin functions */
298 double pow_dd(doublereal *, doublereal *);
299 integer s_cmp(char *, char *, ftnlen, ftnlen);
300 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
301 double d_imag(doublecomplex *);
302 void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
303
304 /* Local variables */
305 static integer j, k, ih, jj, iq, np;
306 static doublecomplex vl[1];
307 static integer wr, ibd, ldh, ldq;
308 static doublereal sep;
309 static integer irz, mode;
310 static doublereal eps23;
311 static integer ierr;
312 static doublecomplex temp;
313 static integer iwev;
314 static char type__[6];
315 static integer ritz, iheig, ihbds;
316 static doublereal conds;
317 static logical reord;
318 extern /* Subroutine */ int zscal_(integer *, doublecomplex *,
319 doublecomplex *, integer *);
320 static integer nconv;
321 extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *,
322 doublecomplex *, integer *, doublecomplex *, integer *);
323 static doublereal rtemp;
324 static doublecomplex rnorm;
325 extern /* Subroutine */ int zgeru_(integer *, integer *, doublecomplex *,
326 doublecomplex *, integer *, doublecomplex *, integer *,
327 doublecomplex *, integer *), zcopy_(integer *, doublecomplex *,
328 integer *, doublecomplex *, integer *), ivout_(integer *, integer
329 *, integer *, integer *, char *, ftnlen), ztrmm_(char *, char *,
330 char *, char *, integer *, integer *, doublecomplex *,
331 doublecomplex *, integer *, doublecomplex *, integer *, ftnlen,
332 ftnlen, ftnlen, ftnlen), zmout_(integer *, integer *, integer *,
333 doublecomplex *, integer *, integer *, char *, ftnlen), zvout_(
334 integer *, integer *, doublecomplex *, integer *, char *, ftnlen);
335 extern doublereal dlapy2_(doublereal *, doublereal *);
336 extern /* Subroutine */ int zgeqr2_(integer *, integer *, doublecomplex *,
337 integer *, doublecomplex *, doublecomplex *, integer *);
338 extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(
339 char *, ftnlen);
340 extern /* Subroutine */ int zunm2r_(char *, char *, integer *, integer *,
341 integer *, doublecomplex *, integer *, doublecomplex *,
342 doublecomplex *, integer *, doublecomplex *, integer *, ftnlen,
343 ftnlen);
344 static integer bounds, invsub, iuptri, msglvl, outncv, numcnv, ishift;
345 extern /* Subroutine */ int zlacpy_(char *, integer *, integer *,
346 doublecomplex *, integer *, doublecomplex *, integer *, ftnlen),
347 zlahqr_(logical *, logical *, integer *, integer *, integer *,
348 doublecomplex *, integer *, doublecomplex *, integer *, integer *,
349 doublecomplex *, integer *, integer *), zngets_(integer *, char *
350 , integer *, integer *, doublecomplex *, doublecomplex *, ftnlen),
351 zlaset_(char *, integer *, integer *, doublecomplex *,
352 doublecomplex *, doublecomplex *, integer *, ftnlen), ztrsen_(
353 char *, char *, logical *, integer *, doublecomplex *, integer *,
354 doublecomplex *, integer *, doublecomplex *, integer *,
355 doublereal *, doublereal *, doublecomplex *, integer *, integer *,
356 ftnlen, ftnlen), ztrevc_(char *, char *, logical *, integer *,
357 doublecomplex *, integer *, doublecomplex *, integer *,
358 doublecomplex *, integer *, integer *, integer *, doublecomplex *,
359 doublereal *, integer *, ftnlen, ftnlen), zdscal_(integer *,
360 doublereal *, doublecomplex *, integer *);
361
362
363 /* %----------------------------------------------------% */
364 /* | Include files for debugging and timing information | */
365 /* %----------------------------------------------------% */
366
367
368 /* \SCCS Information: @(#) */
369 /* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
370
371 /* %---------------------------------% */
372 /* | See debug.doc for documentation | */
373 /* %---------------------------------% */
374
375 /* %------------------% */
376 /* | Scalar Arguments | */
377 /* %------------------% */
378
379 /* %--------------------------------% */
380 /* | See stat.doc for documentation | */
381 /* %--------------------------------% */
382
383 /* \SCCS Information: @(#) */
384 /* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
385
386
387
388 /* %-----------------% */
389 /* | Array Arguments | */
390 /* %-----------------% */
391
392
393 /* %------------% */
394 /* | Parameters | */
395 /* %------------% */
396
397
398 /* %---------------% */
399 /* | Local Scalars | */
400 /* %---------------% */
401
402
403 /* %----------------------% */
404 /* | External Subroutines | */
405 /* %----------------------% */
406
407
408 /* %--------------------% */
409 /* | External Functions | */
410 /* %--------------------% */
411
412
413
414 /* %-----------------------% */
415 /* | Executable Statements | */
416 /* %-----------------------% */
417
418 /* %------------------------% */
419 /* | Set default parameters | */
420 /* %------------------------% */
421
422 /* Parameter adjustments */
423 --workd;
424 --resid;
425 z_dim1 = *ldz;
426 z_offset = 1 + z_dim1;
427 z__ -= z_offset;
428 --d__;
429 --rwork;
430 --workev;
431 --select;
432 v_dim1 = *ldv;
433 v_offset = 1 + v_dim1;
434 v -= v_offset;
435 --iparam;
436 --ipntr;
437 --workl;
438
439 /* Function Body */
440 msglvl = debug_1.mceupd;
441 mode = iparam[7];
442 nconv = iparam[5];
443 *info = 0;
444
445
446 /* %---------------------------------% */
447 /* | Get machine dependent constant. | */
448 /* %---------------------------------% */
449
450 eps23 = dlamch_("Epsilon-Machine", (ftnlen)15);
451 eps23 = pow_dd(&eps23, &c_b5);
452
453 /* %-------------------------------% */
454 /* | Quick return | */
455 /* | Check for incompatible input | */
456 /* %-------------------------------% */
457
458 ierr = 0;
459
460 if (nconv <= 0) {
461 ierr = -14;
462 } else if (*n <= 0) {
463 ierr = -1;
464 } else if (*nev <= 0) {
465 ierr = -2;
466 } else if (*ncv <= *nev || *ncv > *n) {
467 ierr = -3;
468 } else if (s_cmp(which, "LM", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which,
469 "SM", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which, "LR", (ftnlen)2,
470 (ftnlen)2) != 0 && s_cmp(which, "SR", (ftnlen)2, (ftnlen)2) != 0
471 && s_cmp(which, "LI", (ftnlen)2, (ftnlen)2) != 0 && s_cmp(which,
472 "SI", (ftnlen)2, (ftnlen)2) != 0) {
473 ierr = -5;
474 } else if (*(unsigned char *)bmat != 'I' && *(unsigned char *)bmat != 'G')
475 {
476 ierr = -6;
477 } else /* if(complicated condition) */ {
478 /* Computing 2nd power */
479 i__1 = *ncv;
480 if (*lworkl < i__1 * i__1 * 3 + (*ncv << 2)) {
481 ierr = -7;
482 } else if (*(unsigned char *)howmny != 'A' && *(unsigned char *)
483 howmny != 'P' && *(unsigned char *)howmny != 'S' && *rvec) {
484 ierr = -13;
485 } else if (*(unsigned char *)howmny == 'S') {
486 ierr = -12;
487 }
488 }
489
490 if (mode == 1 || mode == 2) {
491 s_copy(type__, "REGULR", (ftnlen)6, (ftnlen)6);
492 } else if (mode == 3) {
493 s_copy(type__, "SHIFTI", (ftnlen)6, (ftnlen)6);
494 } else {
495 ierr = -10;
496 }
497 if (mode == 1 && *(unsigned char *)bmat == 'G') {
498 ierr = -11;
499 }
500
501 /* %------------% */
502 /* | Error Exit | */
503 /* %------------% */
504
505 if (ierr != 0) {
506 *info = ierr;
507 goto L9000;
508 }
509
510 /* %--------------------------------------------------------% */
511 /* | Pointer into WORKL for address of H, RITZ, WORKEV, Q | */
512 /* | etc... and the remaining workspace. | */
513 /* | Also update pointer to be used on output. | */
514 /* | Memory is laid out as follows: | */
515 /* | workl(1:ncv*ncv) := generated Hessenberg matrix | */
516 /* | workl(ncv*ncv+1:ncv*ncv+ncv) := ritz values | */
517 /* | workl(ncv*ncv+ncv+1:ncv*ncv+2*ncv) := error bounds | */
518 /* %--------------------------------------------------------% */
519
520 /* %-----------------------------------------------------------% */
521 /* | The following is used and set by ZNEUPD. | */
522 /* | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := The untransformed | */
523 /* | Ritz values. | */
524 /* | workl(ncv*ncv+3*ncv+1:ncv*ncv+4*ncv) := The untransformed | */
525 /* | error bounds of | */
526 /* | the Ritz values | */
527 /* | workl(ncv*ncv+4*ncv+1:2*ncv*ncv+4*ncv) := Holds the upper | */
528 /* | triangular matrix | */
529 /* | for H. | */
530 /* | workl(2*ncv*ncv+4*ncv+1: 3*ncv*ncv+4*ncv) := Holds the | */
531 /* | associated matrix | */
532 /* | representation of | */
533 /* | the invariant | */
534 /* | subspace for H. | */
535 /* | GRAND total of NCV * ( 3 * NCV + 4 ) locations. | */
536 /* %-----------------------------------------------------------% */
537
538 ih = ipntr[5];
539 ritz = ipntr[6];
540 iq = ipntr[7];
541 bounds = ipntr[8];
542 ldh = *ncv;
543 ldq = *ncv;
544 iheig = bounds + ldh;
545 ihbds = iheig + ldh;
546 iuptri = ihbds + ldh;
547 invsub = iuptri + ldh * *ncv;
548 ipntr[9] = iheig;
549 ipntr[11] = ihbds;
550 ipntr[12] = iuptri;
551 ipntr[13] = invsub;
552 wr = 1;
553 iwev = wr + *ncv;
554
555 /* %-----------------------------------------% */
556 /* | irz points to the Ritz values computed | */
557 /* | by _neigh before exiting _naup2. | */
558 /* | ibd points to the Ritz estimates | */
559 /* | computed by _neigh before exiting | */
560 /* | _naup2. | */
561 /* %-----------------------------------------% */
562
563 irz = ipntr[14] + *ncv * *ncv;
564 ibd = irz + *ncv;
565
566 /* %------------------------------------% */
567 /* | RNORM is B-norm of the RESID(1:N). | */
568 /* %------------------------------------% */
569
570 i__1 = ih + 2;
571 rnorm.r = workl[i__1].r, rnorm.i = workl[i__1].i;
572 i__1 = ih + 2;
573 workl[i__1].r = 0., workl[i__1].i = 0.;
574
575 if (msglvl > 2) {
576 zvout_(&debug_1.logfil, ncv, &workl[irz], &debug_1.ndigit, "_neupd: "
577 "Ritz values passed in from _NAUPD.", (ftnlen)42);
578 zvout_(&debug_1.logfil, ncv, &workl[ibd], &debug_1.ndigit, "_neupd: "
579 "Ritz estimates passed in from _NAUPD.", (ftnlen)45);
580 }
581
582 if (*rvec) {
583
584 reord = FALSE_;
585
586 /* %---------------------------------------------------% */
587 /* | Use the temporary bounds array to store indices | */
588 /* | These will be used to mark the select array later | */
589 /* %---------------------------------------------------% */
590
591 i__1 = *ncv;
592 for (j = 1; j <= i__1; ++j) {
593 i__2 = bounds + j - 1;
594 workl[i__2].r = (doublereal) j, workl[i__2].i = 0.;
595 select[j] = FALSE_;
596 /* L10: */
597 }
598
599 /* %-------------------------------------% */
600 /* | Select the wanted Ritz values. | */
601 /* | Sort the Ritz values so that the | */
602 /* | wanted ones appear at the tailing | */
603 /* | NEV positions of workl(irr) and | */
604 /* | workl(iri). Move the corresponding | */
605 /* | error estimates in workl(ibd) | */
606 /* | accordingly. | */
607 /* %-------------------------------------% */
608
609 np = *ncv - *nev;
610 ishift = 0;
611 zngets_(&ishift, which, nev, &np, &workl[irz], &workl[bounds], (
612 ftnlen)2);
613
614 if (msglvl > 2) {
615 zvout_(&debug_1.logfil, ncv, &workl[irz], &debug_1.ndigit, "_neu"
616 "pd: Ritz values after calling _NGETS.", (ftnlen)41);
617 zvout_(&debug_1.logfil, ncv, &workl[bounds], &debug_1.ndigit,
618 "_neupd: Ritz value indices after calling _NGETS.", (
619 ftnlen)48);
620 }
621
622 /* %-----------------------------------------------------% */
623 /* | Record indices of the converged wanted Ritz values | */
624 /* | Mark the select array for possible reordering | */
625 /* %-----------------------------------------------------% */
626
627 numcnv = 0;
628 i__1 = *ncv;
629 for (j = 1; j <= i__1; ++j) {
630 /* Computing MAX */
631 i__2 = irz + *ncv - j;
632 d__3 = workl[i__2].r;
633 d__4 = d_imag(&workl[irz + *ncv - j]);
634 d__1 = eps23, d__2 = dlapy2_(&d__3, &d__4);
635 rtemp = max(d__1,d__2);
636 i__2 = bounds + *ncv - j;
637 jj = (integer) workl[i__2].r;
638 i__2 = ibd + jj - 1;
639 d__1 = workl[i__2].r;
640 d__2 = d_imag(&workl[ibd + jj - 1]);
641 if (numcnv < nconv && dlapy2_(&d__1, &d__2) <= *tol * rtemp) {
642 select[jj] = TRUE_;
643 ++numcnv;
644 if (jj > *nev) {
645 reord = TRUE_;
646 }
647 }
648 /* L11: */
649 }
650
651 /* %-----------------------------------------------------------% */
652 /* | Check the count (numcnv) of converged Ritz values with | */
653 /* | the number (nconv) reported by dnaupd. If these two | */
654 /* | are different then there has probably been an error | */
655 /* | caused by incorrect passing of the dnaupd data. | */
656 /* %-----------------------------------------------------------% */
657
658 if (msglvl > 2) {
659 ivout_(&debug_1.logfil, &c__1, &numcnv, &debug_1.ndigit, "_neupd"
660 ": Number of specified eigenvalues", (ftnlen)39);
661 ivout_(&debug_1.logfil, &c__1, &nconv, &debug_1.ndigit, "_neupd:"
662 " Number of \"converged\" eigenvalues", (ftnlen)41);
663 }
664
665 if (numcnv != nconv) {
666 *info = -15;
667 goto L9000;
668 }
669
670 /* %-------------------------------------------------------% */
671 /* | Call LAPACK routine zlahqr to compute the Schur form | */
672 /* | of the upper Hessenberg matrix returned by ZNAUPD. | */
673 /* | Make a copy of the upper Hessenberg matrix. | */
674 /* | Initialize the Schur vector matrix Q to the identity. | */
675 /* %-------------------------------------------------------% */
676
677 i__1 = ldh * *ncv;
678 zcopy_(&i__1, &workl[ih], &c__1, &workl[iuptri], &c__1);
679 zlaset_("All", ncv, ncv, &c_b2, &c_b1, &workl[invsub], &ldq, (ftnlen)
680 3);
681 zlahqr_(&c_true, &c_true, ncv, &c__1, ncv, &workl[iuptri], &ldh, &
682 workl[iheig], &c__1, ncv, &workl[invsub], &ldq, &ierr);
683 zcopy_(ncv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds], &c__1);
684
685 if (ierr != 0) {
686 *info = -8;
687 goto L9000;
688 }
689
690 if (msglvl > 1) {
691 zvout_(&debug_1.logfil, ncv, &workl[iheig], &debug_1.ndigit,
692 "_neupd: Eigenvalues of H", (ftnlen)24);
693 zvout_(&debug_1.logfil, ncv, &workl[ihbds], &debug_1.ndigit,
694 "_neupd: Last row of the Schur vector matrix", (ftnlen)43)
695 ;
696 if (msglvl > 3) {
697 zmout_(&debug_1.logfil, ncv, ncv, &workl[iuptri], &ldh, &
698 debug_1.ndigit, "_neupd: The upper triangular matrix "
699 , (ftnlen)36);
700 }
701 }
702
703 if (reord) {
704
705 /* %-----------------------------------------------% */
706 /* | Reorder the computed upper triangular matrix. | */
707 /* %-----------------------------------------------% */
708
709 ztrsen_("None", "V", &select[1], ncv, &workl[iuptri], &ldh, &
710 workl[invsub], &ldq, &workl[iheig], &nconv, &conds, &sep,
711 &workev[1], ncv, &ierr, (ftnlen)4, (ftnlen)1);
712
713 if (ierr == 1) {
714 *info = 1;
715 goto L9000;
716 }
717
718 if (msglvl > 2) {
719 zvout_(&debug_1.logfil, ncv, &workl[iheig], &debug_1.ndigit,
720 "_neupd: Eigenvalues of H--reordered", (ftnlen)35);
721 if (msglvl > 3) {
722 zmout_(&debug_1.logfil, ncv, ncv, &workl[iuptri], &ldq, &
723 debug_1.ndigit, "_neupd: Triangular matrix after"
724 " re-ordering", (ftnlen)43);
725 }
726 }
727
728 }
729
730 /* %---------------------------------------------% */
731 /* | Copy the last row of the Schur basis matrix | */
732 /* | to workl(ihbds). This vector will be used | */
733 /* | to compute the Ritz estimates of converged | */
734 /* | Ritz values. | */
735 /* %---------------------------------------------% */
736
737 zcopy_(ncv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds], &c__1);
738
739 /* %--------------------------------------------% */
740 /* | Place the computed eigenvalues of H into D | */
741 /* | if a spectral transformation was not used. | */
742 /* %--------------------------------------------% */
743
744 if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) == 0) {
745 zcopy_(&nconv, &workl[iheig], &c__1, &d__[1], &c__1);
746 }
747
748 /* %----------------------------------------------------------% */
749 /* | Compute the QR factorization of the matrix representing | */
750 /* | the wanted invariant subspace located in the first NCONV | */
751 /* | columns of workl(invsub,ldq). | */
752 /* %----------------------------------------------------------% */
753
754 zgeqr2_(ncv, &nconv, &workl[invsub], &ldq, &workev[1], &workev[*ncv +
755 1], &ierr);
756
757 /* %--------------------------------------------------------% */
758 /* | * Postmultiply V by Q using zunm2r. | */
759 /* | * Copy the first NCONV columns of VQ into Z. | */
760 /* | * Postmultiply Z by R. | */
761 /* | The N by NCONV matrix Z is now a matrix representation | */
762 /* | of the approximate invariant subspace associated with | */
763 /* | the Ritz values in workl(iheig). The first NCONV | */
764 /* | columns of V are now approximate Schur vectors | */
765 /* | associated with the upper triangular matrix of order | */
766 /* | NCONV in workl(iuptri). | */
767 /* %--------------------------------------------------------% */
768
769 zunm2r_("Right", "Notranspose", n, ncv, &nconv, &workl[invsub], &ldq,
770 &workev[1], &v[v_offset], ldv, &workd[*n + 1], &ierr, (ftnlen)
771 5, (ftnlen)11);
772 zlacpy_("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz, (
773 ftnlen)3);
774
775 i__1 = nconv;
776 for (j = 1; j <= i__1; ++j) {
777
778 /* %---------------------------------------------------% */
779 /* | Perform both a column and row scaling if the | */
780 /* | diagonal element of workl(invsub,ldq) is negative | */
781 /* | I'm lazy and don't take advantage of the upper | */
782 /* | triangular form of workl(iuptri,ldq). | */
783 /* | Note that since Q is orthogonal, R is a diagonal | */
784 /* | matrix consisting of plus or minus ones. | */
785 /* %---------------------------------------------------% */
786
787 i__2 = invsub + (j - 1) * ldq + j - 1;
788 if (workl[i__2].r < 0.) {
789 z__1.r = -1., z__1.i = -0.;
790 zscal_(&nconv, &z__1, &workl[iuptri + j - 1], &ldq);
791 z__1.r = -1., z__1.i = -0.;
792 zscal_(&nconv, &z__1, &workl[iuptri + (j - 1) * ldq], &c__1);
793 }
794
795 /* L20: */
796 }
797
798 if (*(unsigned char *)howmny == 'A') {
799
800 /* %--------------------------------------------% */
801 /* | Compute the NCONV wanted eigenvectors of T | */
802 /* | located in workl(iuptri,ldq). | */
803 /* %--------------------------------------------% */
804
805 i__1 = *ncv;
806 for (j = 1; j <= i__1; ++j) {
807 if (j <= nconv) {
808 select[j] = TRUE_;
809 } else {
810 select[j] = FALSE_;
811 }
812 /* L30: */
813 }
814
815 ztrevc_("Right", "Select", &select[1], ncv, &workl[iuptri], &ldq,
816 vl, &c__1, &workl[invsub], &ldq, ncv, &outncv, &workev[1],
817 &rwork[1], &ierr, (ftnlen)5, (ftnlen)6);
818
819 if (ierr != 0) {
820 *info = -9;
821 goto L9000;
822 }
823
824 /* %------------------------------------------------% */
825 /* | Scale the returning eigenvectors so that their | */
826 /* | Euclidean norms are all one. LAPACK subroutine | */
827 /* | ztrevc returns each eigenvector normalized so | */
828 /* | that the element of largest magnitude has | */
829 /* | magnitude 1. | */
830 /* %------------------------------------------------% */
831
832 i__1 = nconv;
833 for (j = 1; j <= i__1; ++j) {
834 rtemp = dznrm2_(ncv, &workl[invsub + (j - 1) * ldq], &c__1);
835 rtemp = 1. / rtemp;
836 zdscal_(ncv, &rtemp, &workl[invsub + (j - 1) * ldq], &c__1);
837
838 /* %------------------------------------------% */
839 /* | Ritz estimates can be obtained by taking | */
840 /* | the inner product of the last row of the | */
841 /* | Schur basis of H with eigenvectors of T. | */
842 /* | Note that the eigenvector matrix of T is | */
843 /* | upper triangular, thus the length of the | */
844 /* | inner product can be set to j. | */
845 /* %------------------------------------------% */
846
847 i__2 = j;
848 zdotc_(&z__1, &j, &workl[ihbds], &c__1, &workl[invsub + (j -
849 1) * ldq], &c__1);
850 workev[i__2].r = z__1.r, workev[i__2].i = z__1.i;
851 /* L40: */
852 }
853
854 if (msglvl > 2) {
855 zcopy_(&nconv, &workl[invsub + *ncv - 1], &ldq, &workl[ihbds],
856 &c__1);
857 zvout_(&debug_1.logfil, &nconv, &workl[ihbds], &
858 debug_1.ndigit, "_neupd: Last row of the eigenvector"
859 " matrix for T", (ftnlen)48);
860 if (msglvl > 3) {
861 zmout_(&debug_1.logfil, ncv, ncv, &workl[invsub], &ldq, &
862 debug_1.ndigit, "_neupd: The eigenvector matrix "
863 "for T", (ftnlen)36);
864 }
865 }
866
867 /* %---------------------------------------% */
868 /* | Copy Ritz estimates into workl(ihbds) | */
869 /* %---------------------------------------% */
870
871 zcopy_(&nconv, &workev[1], &c__1, &workl[ihbds], &c__1);
872
873 /* %----------------------------------------------% */
874 /* | The eigenvector matrix Q of T is triangular. | */
875 /* | Form Z*Q. | */
876 /* %----------------------------------------------% */
877
878 ztrmm_("Right", "Upper", "No transpose", "Non-unit", n, &nconv, &
879 c_b1, &workl[invsub], &ldq, &z__[z_offset], ldz, (ftnlen)
880 5, (ftnlen)5, (ftnlen)12, (ftnlen)8);
881 }
882
883 } else {
884
885 /* %--------------------------------------------------% */
886 /* | An approximate invariant subspace is not needed. | */
887 /* | Place the Ritz values computed ZNAUPD into D. | */
888 /* %--------------------------------------------------% */
889
890 zcopy_(&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
891 zcopy_(&nconv, &workl[ritz], &c__1, &workl[iheig], &c__1);
892 zcopy_(&nconv, &workl[bounds], &c__1, &workl[ihbds], &c__1);
893
894 }
895
896 /* %------------------------------------------------% */
897 /* | Transform the Ritz values and possibly vectors | */
898 /* | and corresponding error bounds of OP to those | */
899 /* | of A*x = lambda*B*x. | */
900 /* %------------------------------------------------% */
901
902 if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) == 0) {
903
904 if (*rvec) {
905 zscal_(ncv, &rnorm, &workl[ihbds], &c__1);
906 }
907
908 } else {
909
910 /* %---------------------------------------% */
911 /* | A spectral transformation was used. | */
912 /* | * Determine the Ritz estimates of the | */
913 /* | Ritz values in the original system. | */
914 /* %---------------------------------------% */
915
916 if (*rvec) {
917 zscal_(ncv, &rnorm, &workl[ihbds], &c__1);
918 }
919
920 i__1 = *ncv;
921 for (k = 1; k <= i__1; ++k) {
922 i__2 = iheig + k - 1;
923 temp.r = workl[i__2].r, temp.i = workl[i__2].i;
924 i__2 = ihbds + k - 1;
925 z_div(&z__2, &workl[ihbds + k - 1], &temp);
926 z_div(&z__1, &z__2, &temp);
927 workl[i__2].r = z__1.r, workl[i__2].i = z__1.i;
928 /* L50: */
929 }
930
931 }
932
933 /* %-----------------------------------------------------------% */
934 /* | * Transform the Ritz values back to the original system. | */
935 /* | For TYPE = 'SHIFTI' the transformation is | */
936 /* | lambda = 1/theta + sigma | */
937 /* | NOTES: | */
938 /* | *The Ritz vectors are not affected by the transformation. | */
939 /* %-----------------------------------------------------------% */
940
941 if (s_cmp(type__, "SHIFTI", (ftnlen)6, (ftnlen)6) == 0) {
942 i__1 = nconv;
943 for (k = 1; k <= i__1; ++k) {
944 i__2 = k;
945 z_div(&z__2, &c_b1, &workl[iheig + k - 1]);
946 z__1.r = z__2.r + sigma->r, z__1.i = z__2.i + sigma->i;
947 d__[i__2].r = z__1.r, d__[i__2].i = z__1.i;
948 /* L60: */
949 }
950 }
951
952 if (s_cmp(type__, "REGULR", (ftnlen)6, (ftnlen)6) != 0 && msglvl > 1) {
953 zvout_(&debug_1.logfil, &nconv, &d__[1], &debug_1.ndigit, "_neupd: U"
954 "ntransformed Ritz values.", (ftnlen)34);
955 zvout_(&debug_1.logfil, &nconv, &workl[ihbds], &debug_1.ndigit, "_ne"
956 "upd: Ritz estimates of the untransformed Ritz values.", (
957 ftnlen)56);
958 } else if (msglvl > 1) {
959 zvout_(&debug_1.logfil, &nconv, &d__[1], &debug_1.ndigit, "_neupd: C"
960 "onverged Ritz values.", (ftnlen)30);
961 zvout_(&debug_1.logfil, &nconv, &workl[ihbds], &debug_1.ndigit, "_ne"
962 "upd: Associated Ritz estimates.", (ftnlen)34);
963 }
964
965 /* %-------------------------------------------------% */
966 /* | Eigenvector Purification step. Formally perform | */
967 /* | one of inverse subspace iteration. Only used | */
968 /* | for MODE = 3. See reference 3. | */
969 /* %-------------------------------------------------% */
970
971 if (*rvec && *(unsigned char *)howmny == 'A' && s_cmp(type__, "SHIFTI", (
972 ftnlen)6, (ftnlen)6) == 0) {
973
974 /* %------------------------------------------------% */
975 /* | Purify the computed Ritz vectors by adding a | */
976 /* | little bit of the residual vector: | */
977 /* | T | */
978 /* | resid(:)*( e s ) / theta | */
979 /* | NCV | */
980 /* | where H s = s theta. | */
981 /* %------------------------------------------------% */
982
983 i__1 = nconv;
984 for (j = 1; j <= i__1; ++j) {
985 i__2 = iheig + j - 1;
986 if (workl[i__2].r != 0. || workl[i__2].i != 0.) {
987 i__2 = j;
988 z_div(&z__1, &workl[invsub + (j - 1) * ldq + *ncv - 1], &
989 workl[iheig + j - 1]);
990 workev[i__2].r = z__1.r, workev[i__2].i = z__1.i;
991 }
992 /* L100: */
993 }
994 /* %---------------------------------------% */
995 /* | Perform a rank one update to Z and | */
996 /* | purify all the Ritz vectors together. | */
997 /* %---------------------------------------% */
998
999 zgeru_(n, &nconv, &c_b1, &resid[1], &c__1, &workev[1], &c__1, &z__[
1000 z_offset], ldz);
1001
1002 }
1003
1004 L9000:
1005
1006 return 0;
1007
1008 /* %---------------% */
1009 /* | End of zneupd| */
1010 /* %---------------% */
1011
1012 } /* zneupd_ */
1013
1014