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