1 /* ./src_f77/zsteqr.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 doublecomplex c_b1 = {0.,0.};
11 static doublecomplex c_b2 = {1.,0.};
12 static integer c__0 = 0;
13 static integer c__1 = 1;
14 static integer c__2 = 2;
15 static doublereal c_b41 = 1.;
16
zsteqr_(char * compz,integer * n,doublereal * d__,doublereal * e,doublecomplex * z__,integer * ldz,doublereal * work,integer * info,ftnlen compz_len)17 /* Subroutine */ int zsteqr_(char *compz, integer *n, doublereal *d__,
18 doublereal *e, doublecomplex *z__, integer *ldz, doublereal *work,
19 integer *info, ftnlen compz_len)
20 {
21 /* System generated locals */
22 integer z_dim1, z_offset, i__1, i__2;
23 doublereal d__1, d__2;
24
25 /* Builtin functions */
26 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
27
28 /* Local variables */
29 static doublereal b, c__, f, g;
30 static integer i__, j, k, l, m;
31 static doublereal p, r__, s;
32 static integer l1, ii, mm, lm1, mm1, nm1;
33 static doublereal rt1, rt2, eps;
34 static integer lsv;
35 static doublereal tst, eps2;
36 static integer lend, jtot;
37 extern /* Subroutine */ int dlae2_(doublereal *, doublereal *, doublereal
38 *, doublereal *, doublereal *);
39 extern logical lsame_(char *, char *, ftnlen, ftnlen);
40 static doublereal anorm;
41 extern /* Subroutine */ int zlasr_(char *, char *, char *, integer *,
42 integer *, doublereal *, doublereal *, doublecomplex *, integer *,
43 ftnlen, ftnlen, ftnlen), zswap_(integer *, doublecomplex *,
44 integer *, doublecomplex *, integer *), dlaev2_(doublereal *,
45 doublereal *, doublereal *, doublereal *, doublereal *,
46 doublereal *, doublereal *);
47 static integer lendm1, lendp1;
48 extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *,
49 ftnlen);
50 static integer iscale;
51 extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
52 doublereal *, doublereal *, integer *, integer *, doublereal *,
53 integer *, integer *, ftnlen);
54 static doublereal safmin;
55 extern /* Subroutine */ int dlartg_(doublereal *, doublereal *,
56 doublereal *, doublereal *, doublereal *);
57 static doublereal safmax;
58 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
59 extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *,
60 ftnlen);
61 extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *,
62 integer *, ftnlen);
63 static integer lendsv;
64 static doublereal ssfmin;
65 static integer nmaxit, icompz;
66 static doublereal ssfmax;
67 extern /* Subroutine */ int zlaset_(char *, integer *, integer *,
68 doublecomplex *, doublecomplex *, doublecomplex *, integer *,
69 ftnlen);
70
71
72 /* -- LAPACK routine (version 3.0) -- */
73 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
74 /* Courant Institute, Argonne National Lab, and Rice University */
75 /* September 30, 1994 */
76
77 /* .. Scalar Arguments .. */
78 /* .. */
79 /* .. Array Arguments .. */
80 /* .. */
81
82 /* Purpose */
83 /* ======= */
84
85 /* ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a */
86 /* symmetric tridiagonal matrix using the implicit QL or QR method. */
87 /* The eigenvectors of a full or band complex Hermitian matrix can also */
88 /* be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this */
89 /* matrix to tridiagonal form. */
90
91 /* Arguments */
92 /* ========= */
93
94 /* COMPZ (input) CHARACTER*1 */
95 /* = 'N': Compute eigenvalues only. */
96 /* = 'V': Compute eigenvalues and eigenvectors of the original */
97 /* Hermitian matrix. On entry, Z must contain the */
98 /* unitary matrix used to reduce the original matrix */
99 /* to tridiagonal form. */
100 /* = 'I': Compute eigenvalues and eigenvectors of the */
101 /* tridiagonal matrix. Z is initialized to the identity */
102 /* matrix. */
103
104 /* N (input) INTEGER */
105 /* The order of the matrix. N >= 0. */
106
107 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
108 /* On entry, the diagonal elements of the tridiagonal matrix. */
109 /* On exit, if INFO = 0, the eigenvalues in ascending order. */
110
111 /* E (input/output) DOUBLE PRECISION array, dimension (N-1) */
112 /* On entry, the (n-1) subdiagonal elements of the tridiagonal */
113 /* matrix. */
114 /* On exit, E has been destroyed. */
115
116 /* Z (input/output) COMPLEX*16 array, dimension (LDZ, N) */
117 /* On entry, if COMPZ = 'V', then Z contains the unitary */
118 /* matrix used in the reduction to tridiagonal form. */
119 /* On exit, if INFO = 0, then if COMPZ = 'V', Z contains the */
120 /* orthonormal eigenvectors of the original Hermitian matrix, */
121 /* and if COMPZ = 'I', Z contains the orthonormal eigenvectors */
122 /* of the symmetric tridiagonal matrix. */
123 /* If COMPZ = 'N', then Z is not referenced. */
124
125 /* LDZ (input) INTEGER */
126 /* The leading dimension of the array Z. LDZ >= 1, and if */
127 /* eigenvectors are desired, then LDZ >= max(1,N). */
128
129 /* WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) */
130 /* If COMPZ = 'N', then WORK is not referenced. */
131
132 /* INFO (output) INTEGER */
133 /* = 0: successful exit */
134 /* < 0: if INFO = -i, the i-th argument had an illegal value */
135 /* > 0: the algorithm has failed to find all the eigenvalues in */
136 /* a total of 30*N iterations; if INFO = i, then i */
137 /* elements of E have not converged to zero; on exit, D */
138 /* and E contain the elements of a symmetric tridiagonal */
139 /* matrix which is unitarily similar to the original */
140 /* matrix. */
141
142 /* ===================================================================== */
143
144 /* .. Parameters .. */
145 /* .. */
146 /* .. Local Scalars .. */
147 /* .. */
148 /* .. External Functions .. */
149 /* .. */
150 /* .. External Subroutines .. */
151 /* .. */
152 /* .. Intrinsic Functions .. */
153 /* .. */
154 /* .. Executable Statements .. */
155
156 /* Test the input parameters. */
157
158 /* Parameter adjustments */
159 --d__;
160 --e;
161 z_dim1 = *ldz;
162 z_offset = 1 + z_dim1;
163 z__ -= z_offset;
164 --work;
165
166 /* Function Body */
167 *info = 0;
168
169 if (lsame_(compz, "N", (ftnlen)1, (ftnlen)1)) {
170 icompz = 0;
171 } else if (lsame_(compz, "V", (ftnlen)1, (ftnlen)1)) {
172 icompz = 1;
173 } else if (lsame_(compz, "I", (ftnlen)1, (ftnlen)1)) {
174 icompz = 2;
175 } else {
176 icompz = -1;
177 }
178 if (icompz < 0) {
179 *info = -1;
180 } else if (*n < 0) {
181 *info = -2;
182 } else if (*ldz < 1 || icompz > 0 && *ldz < max(1,*n)) {
183 *info = -6;
184 }
185 if (*info != 0) {
186 i__1 = -(*info);
187 xerbla_("ZSTEQR", &i__1, (ftnlen)6);
188 return 0;
189 }
190
191 /* Quick return if possible */
192
193 if (*n == 0) {
194 return 0;
195 }
196
197 if (*n == 1) {
198 if (icompz == 2) {
199 i__1 = z_dim1 + 1;
200 z__[i__1].r = 1., z__[i__1].i = 0.;
201 }
202 return 0;
203 }
204
205 /* Determine the unit roundoff and over/underflow thresholds. */
206
207 eps = dlamch_("E", (ftnlen)1);
208 /* Computing 2nd power */
209 d__1 = eps;
210 eps2 = d__1 * d__1;
211 safmin = dlamch_("S", (ftnlen)1);
212 safmax = 1. / safmin;
213 ssfmax = sqrt(safmax) / 3.;
214 ssfmin = sqrt(safmin) / eps2;
215
216 /* Compute the eigenvalues and eigenvectors of the tridiagonal */
217 /* matrix. */
218
219 if (icompz == 2) {
220 zlaset_("Full", n, n, &c_b1, &c_b2, &z__[z_offset], ldz, (ftnlen)4);
221 }
222
223 nmaxit = *n * 30;
224 jtot = 0;
225
226 /* Determine where the matrix splits and choose QL or QR iteration */
227 /* for each block, according to whether top or bottom diagonal */
228 /* element is smaller. */
229
230 l1 = 1;
231 nm1 = *n - 1;
232
233 L10:
234 if (l1 > *n) {
235 goto L160;
236 }
237 if (l1 > 1) {
238 e[l1 - 1] = 0.;
239 }
240 if (l1 <= nm1) {
241 i__1 = nm1;
242 for (m = l1; m <= i__1; ++m) {
243 tst = (d__1 = e[m], abs(d__1));
244 if (tst == 0.) {
245 goto L30;
246 }
247 if (tst <= sqrt((d__1 = d__[m], abs(d__1))) * sqrt((d__2 = d__[m
248 + 1], abs(d__2))) * eps) {
249 e[m] = 0.;
250 goto L30;
251 }
252 /* L20: */
253 }
254 }
255 m = *n;
256
257 L30:
258 l = l1;
259 lsv = l;
260 lend = m;
261 lendsv = lend;
262 l1 = m + 1;
263 if (lend == l) {
264 goto L10;
265 }
266
267 /* Scale submatrix in rows and columns L to LEND */
268
269 i__1 = lend - l + 1;
270 anorm = dlanst_("I", &i__1, &d__[l], &e[l], (ftnlen)1);
271 iscale = 0;
272 if (anorm == 0.) {
273 goto L10;
274 }
275 if (anorm > ssfmax) {
276 iscale = 1;
277 i__1 = lend - l + 1;
278 dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
279 info, (ftnlen)1);
280 i__1 = lend - l;
281 dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
282 info, (ftnlen)1);
283 } else if (anorm < ssfmin) {
284 iscale = 2;
285 i__1 = lend - l + 1;
286 dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
287 info, (ftnlen)1);
288 i__1 = lend - l;
289 dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
290 info, (ftnlen)1);
291 }
292
293 /* Choose between QL and QR iteration */
294
295 if ((d__1 = d__[lend], abs(d__1)) < (d__2 = d__[l], abs(d__2))) {
296 lend = lsv;
297 l = lendsv;
298 }
299
300 if (lend > l) {
301
302 /* QL Iteration */
303
304 /* Look for small subdiagonal element. */
305
306 L40:
307 if (l != lend) {
308 lendm1 = lend - 1;
309 i__1 = lendm1;
310 for (m = l; m <= i__1; ++m) {
311 /* Computing 2nd power */
312 d__2 = (d__1 = e[m], abs(d__1));
313 tst = d__2 * d__2;
314 if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m
315 + 1], abs(d__2)) + safmin) {
316 goto L60;
317 }
318 /* L50: */
319 }
320 }
321
322 m = lend;
323
324 L60:
325 if (m < lend) {
326 e[m] = 0.;
327 }
328 p = d__[l];
329 if (m == l) {
330 goto L80;
331 }
332
333 /* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 */
334 /* to compute its eigensystem. */
335
336 if (m == l + 1) {
337 if (icompz > 0) {
338 dlaev2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
339 work[l] = c__;
340 work[*n - 1 + l] = s;
341 zlasr_("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
342 z__[l * z_dim1 + 1], ldz, (ftnlen)1, (ftnlen)1, (
343 ftnlen)1);
344 } else {
345 dlae2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
346 }
347 d__[l] = rt1;
348 d__[l + 1] = rt2;
349 e[l] = 0.;
350 l += 2;
351 if (l <= lend) {
352 goto L40;
353 }
354 goto L140;
355 }
356
357 if (jtot == nmaxit) {
358 goto L140;
359 }
360 ++jtot;
361
362 /* Form shift. */
363
364 g = (d__[l + 1] - p) / (e[l] * 2.);
365 r__ = dlapy2_(&g, &c_b41);
366 g = d__[m] - p + e[l] / (g + d_sign(&r__, &g));
367
368 s = 1.;
369 c__ = 1.;
370 p = 0.;
371
372 /* Inner loop */
373
374 mm1 = m - 1;
375 i__1 = l;
376 for (i__ = mm1; i__ >= i__1; --i__) {
377 f = s * e[i__];
378 b = c__ * e[i__];
379 dlartg_(&g, &f, &c__, &s, &r__);
380 if (i__ != m - 1) {
381 e[i__ + 1] = r__;
382 }
383 g = d__[i__ + 1] - p;
384 r__ = (d__[i__] - g) * s + c__ * 2. * b;
385 p = s * r__;
386 d__[i__ + 1] = g + p;
387 g = c__ * r__ - b;
388
389 /* If eigenvectors are desired, then save rotations. */
390
391 if (icompz > 0) {
392 work[i__] = c__;
393 work[*n - 1 + i__] = -s;
394 }
395
396 /* L70: */
397 }
398
399 /* If eigenvectors are desired, then apply saved rotations. */
400
401 if (icompz > 0) {
402 mm = m - l + 1;
403 zlasr_("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l
404 * z_dim1 + 1], ldz, (ftnlen)1, (ftnlen)1, (ftnlen)1);
405 }
406
407 d__[l] -= p;
408 e[l] = g;
409 goto L40;
410
411 /* Eigenvalue found. */
412
413 L80:
414 d__[l] = p;
415
416 ++l;
417 if (l <= lend) {
418 goto L40;
419 }
420 goto L140;
421
422 } else {
423
424 /* QR Iteration */
425
426 /* Look for small superdiagonal element. */
427
428 L90:
429 if (l != lend) {
430 lendp1 = lend + 1;
431 i__1 = lendp1;
432 for (m = l; m >= i__1; --m) {
433 /* Computing 2nd power */
434 d__2 = (d__1 = e[m - 1], abs(d__1));
435 tst = d__2 * d__2;
436 if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m
437 - 1], abs(d__2)) + safmin) {
438 goto L110;
439 }
440 /* L100: */
441 }
442 }
443
444 m = lend;
445
446 L110:
447 if (m > lend) {
448 e[m - 1] = 0.;
449 }
450 p = d__[l];
451 if (m == l) {
452 goto L130;
453 }
454
455 /* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 */
456 /* to compute its eigensystem. */
457
458 if (m == l - 1) {
459 if (icompz > 0) {
460 dlaev2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
461 ;
462 work[m] = c__;
463 work[*n - 1 + m] = s;
464 zlasr_("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
465 z__[(l - 1) * z_dim1 + 1], ldz, (ftnlen)1, (ftnlen)1,
466 (ftnlen)1);
467 } else {
468 dlae2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
469 }
470 d__[l - 1] = rt1;
471 d__[l] = rt2;
472 e[l - 1] = 0.;
473 l += -2;
474 if (l >= lend) {
475 goto L90;
476 }
477 goto L140;
478 }
479
480 if (jtot == nmaxit) {
481 goto L140;
482 }
483 ++jtot;
484
485 /* Form shift. */
486
487 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
488 r__ = dlapy2_(&g, &c_b41);
489 g = d__[m] - p + e[l - 1] / (g + d_sign(&r__, &g));
490
491 s = 1.;
492 c__ = 1.;
493 p = 0.;
494
495 /* Inner loop */
496
497 lm1 = l - 1;
498 i__1 = lm1;
499 for (i__ = m; i__ <= i__1; ++i__) {
500 f = s * e[i__];
501 b = c__ * e[i__];
502 dlartg_(&g, &f, &c__, &s, &r__);
503 if (i__ != m) {
504 e[i__ - 1] = r__;
505 }
506 g = d__[i__] - p;
507 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
508 p = s * r__;
509 d__[i__] = g + p;
510 g = c__ * r__ - b;
511
512 /* If eigenvectors are desired, then save rotations. */
513
514 if (icompz > 0) {
515 work[i__] = c__;
516 work[*n - 1 + i__] = s;
517 }
518
519 /* L120: */
520 }
521
522 /* If eigenvectors are desired, then apply saved rotations. */
523
524 if (icompz > 0) {
525 mm = l - m + 1;
526 zlasr_("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m
527 * z_dim1 + 1], ldz, (ftnlen)1, (ftnlen)1, (ftnlen)1);
528 }
529
530 d__[l] -= p;
531 e[lm1] = g;
532 goto L90;
533
534 /* Eigenvalue found. */
535
536 L130:
537 d__[l] = p;
538
539 --l;
540 if (l >= lend) {
541 goto L90;
542 }
543 goto L140;
544
545 }
546
547 /* Undo scaling if necessary */
548
549 L140:
550 if (iscale == 1) {
551 i__1 = lendsv - lsv + 1;
552 dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
553 n, info, (ftnlen)1);
554 i__1 = lendsv - lsv;
555 dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
556 info, (ftnlen)1);
557 } else if (iscale == 2) {
558 i__1 = lendsv - lsv + 1;
559 dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
560 n, info, (ftnlen)1);
561 i__1 = lendsv - lsv;
562 dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
563 info, (ftnlen)1);
564 }
565
566 /* Check for no convergence to an eigenvalue after a total */
567 /* of N*MAXIT iterations. */
568
569 if (jtot == nmaxit) {
570 i__1 = *n - 1;
571 for (i__ = 1; i__ <= i__1; ++i__) {
572 if (e[i__] != 0.) {
573 ++(*info);
574 }
575 /* L150: */
576 }
577 return 0;
578 }
579 goto L10;
580
581 /* Order eigenvalues and eigenvectors. */
582
583 L160:
584 if (icompz == 0) {
585
586 /* Use Quick Sort */
587
588 dlasrt_("I", n, &d__[1], info, (ftnlen)1);
589
590 } else {
591
592 /* Use Selection Sort to minimize swaps of eigenvectors */
593
594 i__1 = *n;
595 for (ii = 2; ii <= i__1; ++ii) {
596 i__ = ii - 1;
597 k = i__;
598 p = d__[i__];
599 i__2 = *n;
600 for (j = ii; j <= i__2; ++j) {
601 if (d__[j] < p) {
602 k = j;
603 p = d__[j];
604 }
605 /* L170: */
606 }
607 if (k != i__) {
608 d__[k] = d__[i__];
609 d__[i__] = p;
610 zswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1],
611 &c__1);
612 }
613 /* L180: */
614 }
615 }
616 return 0;
617
618 /* End of ZSTEQR */
619
620 } /* zsteqr_ */
621
622