1 /* ./src_f77/zhsein.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 /* Table of constant values */
9
10 static logical c_false = FALSE_;
11 static logical c_true = TRUE_;
12
zhsein_(char * side,char * eigsrc,char * initv,logical * select,integer * n,doublecomplex * h__,integer * ldh,doublecomplex * w,doublecomplex * vl,integer * ldvl,doublecomplex * vr,integer * ldvr,integer * mm,integer * m,doublecomplex * work,doublereal * rwork,integer * ifaill,integer * ifailr,integer * info,ftnlen side_len,ftnlen eigsrc_len,ftnlen initv_len)13 /* Subroutine */ int zhsein_(char *side, char *eigsrc, char *initv, logical *
14 select, integer *n, doublecomplex *h__, integer *ldh, doublecomplex *
15 w, doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr,
16 integer *mm, integer *m, doublecomplex *work, doublereal *rwork,
17 integer *ifaill, integer *ifailr, integer *info, ftnlen side_len,
18 ftnlen eigsrc_len, ftnlen initv_len)
19 {
20 /* System generated locals */
21 integer h_dim1, h_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1,
22 i__2, i__3;
23 doublereal d__1, d__2;
24 doublecomplex z__1, z__2;
25
26 /* Builtin functions */
27 double d_imag(doublecomplex *);
28
29 /* Local variables */
30 static integer i__, k, kl, kr, ks;
31 static doublecomplex wk;
32 static integer kln;
33 static doublereal ulp, eps3, unfl;
34 extern logical lsame_(char *, char *, ftnlen, ftnlen);
35 static integer iinfo;
36 static logical leftv, bothv;
37 static doublereal hnorm;
38 extern doublereal dlamch_(char *, ftnlen);
39 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), zlaein_(
40 logical *, logical *, integer *, doublecomplex *, integer *,
41 doublecomplex *, doublecomplex *, doublecomplex *, integer *,
42 doublereal *, doublereal *, doublereal *, integer *);
43 extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *,
44 doublereal *, ftnlen);
45 static logical noinit;
46 static integer ldwork;
47 static logical rightv, fromqr;
48 static doublereal smlnum;
49
50
51 /* -- LAPACK routine (version 3.0) -- */
52 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
53 /* Courant Institute, Argonne National Lab, and Rice University */
54 /* September 30, 1994 */
55
56 /* .. Scalar Arguments .. */
57 /* .. */
58 /* .. Array Arguments .. */
59 /* .. */
60
61 /* Purpose */
62 /* ======= */
63
64 /* ZHSEIN uses inverse iteration to find specified right and/or left */
65 /* eigenvectors of a complex upper Hessenberg matrix H. */
66
67 /* The right eigenvector x and the left eigenvector y of the matrix H */
68 /* corresponding to an eigenvalue w are defined by: */
69
70 /* H * x = w * x, y**h * H = w * y**h */
71
72 /* where y**h denotes the conjugate transpose of the vector y. */
73
74 /* Arguments */
75 /* ========= */
76
77 /* SIDE (input) CHARACTER*1 */
78 /* = 'R': compute right eigenvectors only; */
79 /* = 'L': compute left eigenvectors only; */
80 /* = 'B': compute both right and left eigenvectors. */
81
82 /* EIGSRC (input) CHARACTER*1 */
83 /* Specifies the source of eigenvalues supplied in W: */
84 /* = 'Q': the eigenvalues were found using ZHSEQR; thus, if */
85 /* H has zero subdiagonal elements, and so is */
86 /* block-triangular, then the j-th eigenvalue can be */
87 /* assumed to be an eigenvalue of the block containing */
88 /* the j-th row/column. This property allows ZHSEIN to */
89 /* perform inverse iteration on just one diagonal block. */
90 /* = 'N': no assumptions are made on the correspondence */
91 /* between eigenvalues and diagonal blocks. In this */
92 /* case, ZHSEIN must always perform inverse iteration */
93 /* using the whole matrix H. */
94
95 /* INITV (input) CHARACTER*1 */
96 /* = 'N': no initial vectors are supplied; */
97 /* = 'U': user-supplied initial vectors are stored in the arrays */
98 /* VL and/or VR. */
99
100 /* SELECT (input) LOGICAL array, dimension (N) */
101 /* Specifies the eigenvectors to be computed. To select the */
102 /* eigenvector corresponding to the eigenvalue W(j), */
103 /* SELECT(j) must be set to .TRUE.. */
104
105 /* N (input) INTEGER */
106 /* The order of the matrix H. N >= 0. */
107
108 /* H (input) COMPLEX*16 array, dimension (LDH,N) */
109 /* The upper Hessenberg matrix H. */
110
111 /* LDH (input) INTEGER */
112 /* The leading dimension of the array H. LDH >= max(1,N). */
113
114 /* W (input/output) COMPLEX*16 array, dimension (N) */
115 /* On entry, the eigenvalues of H. */
116 /* On exit, the real parts of W may have been altered since */
117 /* close eigenvalues are perturbed slightly in searching for */
118 /* independent eigenvectors. */
119
120 /* VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) */
121 /* On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must */
122 /* contain starting vectors for the inverse iteration for the */
123 /* left eigenvectors; the starting vector for each eigenvector */
124 /* must be in the same column in which the eigenvector will be */
125 /* stored. */
126 /* On exit, if SIDE = 'L' or 'B', the left eigenvectors */
127 /* specified by SELECT will be stored consecutively in the */
128 /* columns of VL, in the same order as their eigenvalues. */
129 /* If SIDE = 'R', VL is not referenced. */
130
131 /* LDVL (input) INTEGER */
132 /* The leading dimension of the array VL. */
133 /* LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. */
134
135 /* VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) */
136 /* On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must */
137 /* contain starting vectors for the inverse iteration for the */
138 /* right eigenvectors; the starting vector for each eigenvector */
139 /* must be in the same column in which the eigenvector will be */
140 /* stored. */
141 /* On exit, if SIDE = 'R' or 'B', the right eigenvectors */
142 /* specified by SELECT will be stored consecutively in the */
143 /* columns of VR, in the same order as their eigenvalues. */
144 /* If SIDE = 'L', VR is not referenced. */
145
146 /* LDVR (input) INTEGER */
147 /* The leading dimension of the array VR. */
148 /* LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. */
149
150 /* MM (input) INTEGER */
151 /* The number of columns in the arrays VL and/or VR. MM >= M. */
152
153 /* M (output) INTEGER */
154 /* The number of columns in the arrays VL and/or VR required to */
155 /* store the eigenvectors (= the number of .TRUE. elements in */
156 /* SELECT). */
157
158 /* WORK (workspace) COMPLEX*16 array, dimension (N*N) */
159
160 /* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */
161
162 /* IFAILL (output) INTEGER array, dimension (MM) */
163 /* If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left */
164 /* eigenvector in the i-th column of VL (corresponding to the */
165 /* eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the */
166 /* eigenvector converged satisfactorily. */
167 /* If SIDE = 'R', IFAILL is not referenced. */
168
169 /* IFAILR (output) INTEGER array, dimension (MM) */
170 /* If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right */
171 /* eigenvector in the i-th column of VR (corresponding to the */
172 /* eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the */
173 /* eigenvector converged satisfactorily. */
174 /* If SIDE = 'L', IFAILR is not referenced. */
175
176 /* INFO (output) INTEGER */
177 /* = 0: successful exit */
178 /* < 0: if INFO = -i, the i-th argument had an illegal value */
179 /* > 0: if INFO = i, i is the number of eigenvectors which */
180 /* failed to converge; see IFAILL and IFAILR for further */
181 /* details. */
182
183 /* Further Details */
184 /* =============== */
185
186 /* Each eigenvector is normalized so that the element of largest */
187 /* magnitude has magnitude 1; here the magnitude of a complex number */
188 /* (x,y) is taken to be |x|+|y|. */
189
190 /* ===================================================================== */
191
192 /* .. Parameters .. */
193 /* .. */
194 /* .. Local Scalars .. */
195 /* .. */
196 /* .. External Functions .. */
197 /* .. */
198 /* .. External Subroutines .. */
199 /* .. */
200 /* .. Intrinsic Functions .. */
201 /* .. */
202 /* .. Statement Functions .. */
203 /* .. */
204 /* .. Statement Function definitions .. */
205 /* .. */
206 /* .. Executable Statements .. */
207
208 /* Decode and test the input parameters. */
209
210 /* Parameter adjustments */
211 --select;
212 h_dim1 = *ldh;
213 h_offset = 1 + h_dim1;
214 h__ -= h_offset;
215 --w;
216 vl_dim1 = *ldvl;
217 vl_offset = 1 + vl_dim1;
218 vl -= vl_offset;
219 vr_dim1 = *ldvr;
220 vr_offset = 1 + vr_dim1;
221 vr -= vr_offset;
222 --work;
223 --rwork;
224 --ifaill;
225 --ifailr;
226
227 /* Function Body */
228 bothv = lsame_(side, "B", (ftnlen)1, (ftnlen)1);
229 rightv = lsame_(side, "R", (ftnlen)1, (ftnlen)1) || bothv;
230 leftv = lsame_(side, "L", (ftnlen)1, (ftnlen)1) || bothv;
231
232 fromqr = lsame_(eigsrc, "Q", (ftnlen)1, (ftnlen)1);
233
234 noinit = lsame_(initv, "N", (ftnlen)1, (ftnlen)1);
235
236 /* Set M to the number of columns required to store the selected */
237 /* eigenvectors. */
238
239 *m = 0;
240 i__1 = *n;
241 for (k = 1; k <= i__1; ++k) {
242 if (select[k]) {
243 ++(*m);
244 }
245 /* L10: */
246 }
247
248 *info = 0;
249 if (! rightv && ! leftv) {
250 *info = -1;
251 } else if (! fromqr && ! lsame_(eigsrc, "N", (ftnlen)1, (ftnlen)1)) {
252 *info = -2;
253 } else if (! noinit && ! lsame_(initv, "U", (ftnlen)1, (ftnlen)1)) {
254 *info = -3;
255 } else if (*n < 0) {
256 *info = -5;
257 } else if (*ldh < max(1,*n)) {
258 *info = -7;
259 } else if (*ldvl < 1 || leftv && *ldvl < *n) {
260 *info = -10;
261 } else if (*ldvr < 1 || rightv && *ldvr < *n) {
262 *info = -12;
263 } else if (*mm < *m) {
264 *info = -13;
265 }
266 if (*info != 0) {
267 i__1 = -(*info);
268 xerbla_("ZHSEIN", &i__1, (ftnlen)6);
269 return 0;
270 }
271
272 /* Quick return if possible. */
273
274 if (*n == 0) {
275 return 0;
276 }
277
278 /* Set machine-dependent constants. */
279
280 unfl = dlamch_("Safe minimum", (ftnlen)12);
281 ulp = dlamch_("Precision", (ftnlen)9);
282 smlnum = unfl * (*n / ulp);
283
284 ldwork = *n;
285
286 kl = 1;
287 kln = 0;
288 if (fromqr) {
289 kr = 0;
290 } else {
291 kr = *n;
292 }
293 ks = 1;
294
295 i__1 = *n;
296 for (k = 1; k <= i__1; ++k) {
297 if (select[k]) {
298
299 /* Compute eigenvector(s) corresponding to W(K). */
300
301 if (fromqr) {
302
303 /* If affiliation of eigenvalues is known, check whether */
304 /* the matrix splits. */
305
306 /* Determine KL and KR such that 1 <= KL <= K <= KR <= N */
307 /* and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or */
308 /* KR = N). */
309
310 /* Then inverse iteration can be performed with the */
311 /* submatrix H(KL:N,KL:N) for a left eigenvector, and with */
312 /* the submatrix H(1:KR,1:KR) for a right eigenvector. */
313
314 i__2 = kl + 1;
315 for (i__ = k; i__ >= i__2; --i__) {
316 i__3 = i__ + (i__ - 1) * h_dim1;
317 if (h__[i__3].r == 0. && h__[i__3].i == 0.) {
318 goto L30;
319 }
320 /* L20: */
321 }
322 L30:
323 kl = i__;
324 if (k > kr) {
325 i__2 = *n - 1;
326 for (i__ = k; i__ <= i__2; ++i__) {
327 i__3 = i__ + 1 + i__ * h_dim1;
328 if (h__[i__3].r == 0. && h__[i__3].i == 0.) {
329 goto L50;
330 }
331 /* L40: */
332 }
333 L50:
334 kr = i__;
335 }
336 }
337
338 if (kl != kln) {
339 kln = kl;
340
341 /* Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it */
342 /* has not ben computed before. */
343
344 i__2 = kr - kl + 1;
345 hnorm = zlanhs_("I", &i__2, &h__[kl + kl * h_dim1], ldh, &
346 rwork[1], (ftnlen)1);
347 if (hnorm > 0.) {
348 eps3 = hnorm * ulp;
349 } else {
350 eps3 = smlnum;
351 }
352 }
353
354 /* Perturb eigenvalue if it is close to any previous */
355 /* selected eigenvalues affiliated to the submatrix */
356 /* H(KL:KR,KL:KR). Close roots are modified by EPS3. */
357
358 i__2 = k;
359 wk.r = w[i__2].r, wk.i = w[i__2].i;
360 L60:
361 i__2 = kl;
362 for (i__ = k - 1; i__ >= i__2; --i__) {
363 i__3 = i__;
364 z__2.r = w[i__3].r - wk.r, z__2.i = w[i__3].i - wk.i;
365 z__1.r = z__2.r, z__1.i = z__2.i;
366 if (select[i__] && (d__1 = z__1.r, abs(d__1)) + (d__2 =
367 d_imag(&z__1), abs(d__2)) < eps3) {
368 z__1.r = wk.r + eps3, z__1.i = wk.i;
369 wk.r = z__1.r, wk.i = z__1.i;
370 goto L60;
371 }
372 /* L70: */
373 }
374 i__2 = k;
375 w[i__2].r = wk.r, w[i__2].i = wk.i;
376
377 if (leftv) {
378
379 /* Compute left eigenvector. */
380
381 i__2 = *n - kl + 1;
382 zlaein_(&c_false, &noinit, &i__2, &h__[kl + kl * h_dim1], ldh,
383 &wk, &vl[kl + ks * vl_dim1], &work[1], &ldwork, &
384 rwork[1], &eps3, &smlnum, &iinfo);
385 if (iinfo > 0) {
386 ++(*info);
387 ifaill[ks] = k;
388 } else {
389 ifaill[ks] = 0;
390 }
391 i__2 = kl - 1;
392 for (i__ = 1; i__ <= i__2; ++i__) {
393 i__3 = i__ + ks * vl_dim1;
394 vl[i__3].r = 0., vl[i__3].i = 0.;
395 /* L80: */
396 }
397 }
398 if (rightv) {
399
400 /* Compute right eigenvector. */
401
402 zlaein_(&c_true, &noinit, &kr, &h__[h_offset], ldh, &wk, &vr[
403 ks * vr_dim1 + 1], &work[1], &ldwork, &rwork[1], &
404 eps3, &smlnum, &iinfo);
405 if (iinfo > 0) {
406 ++(*info);
407 ifailr[ks] = k;
408 } else {
409 ifailr[ks] = 0;
410 }
411 i__2 = *n;
412 for (i__ = kr + 1; i__ <= i__2; ++i__) {
413 i__3 = i__ + ks * vr_dim1;
414 vr[i__3].r = 0., vr[i__3].i = 0.;
415 /* L90: */
416 }
417 }
418 ++ks;
419 }
420 /* L100: */
421 }
422
423 return 0;
424
425 /* End of ZHSEIN */
426
427 } /* zhsein_ */
428
429