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