3 /* Table of constant values */
5 static doublereal c_b9 = 0.;
6 static doublereal c_b10 = 1.;
7 static integer c__0 = 0;
8 static integer c__1 = 1;
9 static integer c__2 = 2;
11 /* Subroutine */ int dsteqr_(char *compz, integer *n, doublereal *d__,
12 doublereal *e, doublereal *z__, integer *ldz, doublereal *work,
15 /* System generated locals */
16 integer z_dim1, z_offset, i__1, i__2;
17 doublereal d__1, d__2;
19 /* Builtin functions */
20 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
23 doublereal b, c__, f, g;
24 integer i__, j, k, l, m;
26 integer l1, ii, mm, lm1, mm1, nm1;
27 doublereal rt1, rt2, eps;
31 extern /* Subroutine */ int dlae2_(doublereal *, doublereal *, doublereal
32 *, doublereal *, doublereal *);
33 extern logical lsame_(char *, char *);
34 extern /* Subroutine */ int dlasr_(char *, char *, char *, integer *,
35 integer *, doublereal *, doublereal *, doublereal *, integer *);
37 extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
38 doublereal *, integer *), dlaev2_(doublereal *, doublereal *,
39 doublereal *, doublereal *, doublereal *, doublereal *,
41 integer lendm1, lendp1;
42 extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);
44 extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
45 doublereal *, doublereal *, integer *, integer *, doublereal *,
46 integer *, integer *), dlaset_(char *, integer *, integer
47 *, doublereal *, doublereal *, doublereal *, integer *);
49 extern /* Subroutine */ int dlartg_(doublereal *, doublereal *,
50 doublereal *, doublereal *, doublereal *);
52 extern /* Subroutine */ int xerbla_(char *, integer *);
53 extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);
54 extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *,
58 integer nmaxit, icompz;
62 /* -- LAPACK routine (version 3.1) -- */
63 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
66 /* .. Scalar Arguments .. */
68 /* .. Array Arguments .. */
74 /* DSTEQR computes all eigenvalues and, optionally, eigenvectors of a */
75 /* symmetric tridiagonal matrix using the implicit QL or QR method. */
76 /* The eigenvectors of a full or band symmetric matrix can also be found */
77 /* if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to */
78 /* tridiagonal form. */
83 /* COMPZ (input) CHARACTER*1 */
84 /* = 'N': Compute eigenvalues only. */
85 /* = 'V': Compute eigenvalues and eigenvectors of the original */
86 /* symmetric matrix. On entry, Z must contain the */
87 /* orthogonal matrix used to reduce the original matrix */
88 /* to tridiagonal form. */
89 /* = 'I': Compute eigenvalues and eigenvectors of the */
90 /* tridiagonal matrix. Z is initialized to the identity */
93 /* N (input) INTEGER */
94 /* The order of the matrix. N >= 0. */
96 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
97 /* On entry, the diagonal elements of the tridiagonal matrix. */
98 /* On exit, if INFO = 0, the eigenvalues in ascending order. */
100 /* E (input/output) DOUBLE PRECISION array, dimension (N-1) */
101 /* On entry, the (n-1) subdiagonal elements of the tridiagonal */
103 /* On exit, E has been destroyed. */
105 /* Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) */
106 /* On entry, if COMPZ = 'V', then Z contains the orthogonal */
107 /* matrix used in the reduction to tridiagonal form. */
108 /* On exit, if INFO = 0, then if COMPZ = 'V', Z contains the */
109 /* orthonormal eigenvectors of the original symmetric matrix, */
110 /* and if COMPZ = 'I', Z contains the orthonormal eigenvectors */
111 /* of the symmetric tridiagonal matrix. */
112 /* If COMPZ = 'N', then Z is not referenced. */
114 /* LDZ (input) INTEGER */
115 /* The leading dimension of the array Z. LDZ >= 1, and if */
116 /* eigenvectors are desired, then LDZ >= max(1,N). */
118 /* WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) */
119 /* If COMPZ = 'N', then WORK is not referenced. */
121 /* INFO (output) INTEGER */
122 /* = 0: successful exit */
123 /* < 0: if INFO = -i, the i-th argument had an illegal value */
124 /* > 0: the algorithm has failed to find all the eigenvalues in */
125 /* a total of 30*N iterations; if INFO = i, then i */
126 /* elements of E have not converged to zero; on exit, D */
127 /* and E contain the elements of a symmetric tridiagonal */
128 /* matrix which is orthogonally similar to the original */
131 /* ===================================================================== */
133 /* .. Parameters .. */
135 /* .. Local Scalars .. */
137 /* .. External Functions .. */
139 /* .. External Subroutines .. */
141 /* .. Intrinsic Functions .. */
143 /* .. Executable Statements .. */
145 /* Test the input parameters. */
147 /* Parameter adjustments */
151 z_offset = 1 + z_dim1;
158 if (lsame_(compz, "N")) {
160 } else if (lsame_(compz, "V")) {
162 } else if (lsame_(compz, "I")) {
171 } else if (*ldz < 1 || icompz > 0 && *ldz < max(1,*n)) {
176 xerbla_("DSTEQR", &i__1);
180 /* Quick return if possible */
188 z__[z_dim1 + 1] = 1.;
193 /* Determine the unit roundoff and over/underflow thresholds. */
196 /* Computing 2nd power */
199 safmin = dlamch_("S");
200 safmax = 1. / safmin;
201 ssfmax = sqrt(safmax) / 3.;
202 ssfmin = sqrt(safmin) / eps2;
204 /* Compute the eigenvalues and eigenvectors of the tridiagonal */
208 dlaset_("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
214 /* Determine where the matrix splits and choose QL or QR iteration */
215 /* for each block, according to whether top or bottom diagonal */
216 /* element is smaller. */
230 for (m = l1; m <= i__1; ++m) {
231 tst = (d__1 = e[m], abs(d__1));
235 if (tst <= sqrt((d__1 = d__[m], abs(d__1))) * sqrt((d__2 = d__[m
236 + 1], abs(d__2))) * eps) {
255 /* Scale submatrix in rows and columns L to LEND */
258 anorm = dlanst_("I", &i__1, &d__[l], &e[l]);
263 if (anorm > ssfmax) {
266 dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
269 dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
271 } else if (anorm < ssfmin) {
274 dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
277 dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
281 /* Choose between QL and QR iteration */
283 if ((d__1 = d__[lend], abs(d__1)) < (d__2 = d__[l], abs(d__2))) {
292 /* Look for small subdiagonal element. */
298 for (m = l; m <= i__1; ++m) {
299 /* Computing 2nd power */
300 d__2 = (d__1 = e[m], abs(d__1));
302 if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m
303 + 1], abs(d__2)) + safmin) {
321 /* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 */
322 /* to compute its eigensystem. */
326 dlaev2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
328 work[*n - 1 + l] = s;
329 dlasr_("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
330 z__[l * z_dim1 + 1], ldz);
332 dlae2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
344 if (jtot == nmaxit) {
351 g = (d__[l + 1] - p) / (e[l] * 2.);
352 r__ = dlapy2_(&g, &c_b10);
353 g = d__[m] - p + e[l] / (g + d_sign(&r__, &g));
363 for (i__ = mm1; i__ >= i__1; --i__) {
366 dlartg_(&g, &f, &c__, &s, &r__);
370 g = d__[i__ + 1] - p;
371 r__ = (d__[i__] - g) * s + c__ * 2. * b;
373 d__[i__ + 1] = g + p;
376 /* If eigenvectors are desired, then save rotations. */
380 work[*n - 1 + i__] = -s;
386 /* If eigenvectors are desired, then apply saved rotations. */
390 dlasr_("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l
398 /* Eigenvalue found. */
413 /* Look for small superdiagonal element. */
419 for (m = l; m >= i__1; --m) {
420 /* Computing 2nd power */
421 d__2 = (d__1 = e[m - 1], abs(d__1));
423 if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m
424 - 1], abs(d__2)) + safmin) {
442 /* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 */
443 /* to compute its eigensystem. */
447 dlaev2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
450 work[*n - 1 + m] = s;
451 dlasr_("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
452 z__[(l - 1) * z_dim1 + 1], ldz);
454 dlae2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
466 if (jtot == nmaxit) {
473 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
474 r__ = dlapy2_(&g, &c_b10);
475 g = d__[m] - p + e[l - 1] / (g + d_sign(&r__, &g));
485 for (i__ = m; i__ <= i__1; ++i__) {
488 dlartg_(&g, &f, &c__, &s, &r__);
493 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
498 /* If eigenvectors are desired, then save rotations. */
502 work[*n - 1 + i__] = s;
508 /* If eigenvectors are desired, then apply saved rotations. */
512 dlasr_("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m
520 /* Eigenvalue found. */
533 /* Undo scaling if necessary */
537 i__1 = lendsv - lsv + 1;
538 dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
541 dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
543 } else if (iscale == 2) {
544 i__1 = lendsv - lsv + 1;
545 dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
548 dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
552 /* Check for no convergence to an eigenvalue after a total */
553 /* of N*MAXIT iterations. */
559 for (i__ = 1; i__ <= i__1; ++i__) {
567 /* Order eigenvalues and eigenvectors. */
574 dlasrt_("I", n, &d__[1], info);
578 /* Use Selection Sort to minimize swaps of eigenvectors */
581 for (ii = 2; ii <= i__1; ++ii) {
586 for (j = ii; j <= i__2; ++j) {
596 dswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1],