3 /* Table of constant values */
5 static real c_b9 = 0.f;
6 static real c_b10 = 1.f;
7 static integer c__0 = 0;
8 static integer c__1 = 1;
9 static integer c__2 = 2;
11 /* Subroutine */ int ssteqr_(char *compz, integer *n, real *d__, real *e,
12 real *z__, integer *ldz, real *work, integer *info)
14 /* System generated locals */
15 integer z_dim1, z_offset, i__1, i__2;
18 /* Builtin functions */
19 double sqrt(doublereal), r_sign(real *, real *);
23 integer i__, j, k, l, m;
25 integer l1, ii, mm, lm1, mm1, nm1;
30 extern /* Subroutine */ int slae2_(real *, real *, real *, real *, real *)
32 extern logical lsame_(char *, char *);
34 extern /* Subroutine */ int slasr_(char *, char *, char *, integer *,
35 integer *, real *, real *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer *);
36 integer lendm1, lendp1;
37 extern /* Subroutine */ int slaev2_(real *, real *, real *, real *, real *
39 extern doublereal slapy2_(real *, real *);
41 extern doublereal slamch_(char *);
43 extern /* Subroutine */ int xerbla_(char *, integer *);
45 extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
46 real *, integer *, integer *, real *, integer *, integer *);
48 extern /* Subroutine */ int slartg_(real *, real *, real *, real *, real *
49 ), slaset_(char *, integer *, integer *, real *, real *, real *,
52 integer nmaxit, icompz;
54 extern doublereal slanst_(char *, integer *, real *, real *);
55 extern /* Subroutine */ int slasrt_(char *, integer *, real *, integer *);
58 /* -- LAPACK routine (version 3.1) -- */
59 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
62 /* .. Scalar Arguments .. */
64 /* .. Array Arguments .. */
70 /* SSTEQR computes all eigenvalues and, optionally, eigenvectors of a */
71 /* symmetric tridiagonal matrix using the implicit QL or QR method. */
72 /* The eigenvectors of a full or band symmetric matrix can also be found */
73 /* if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to */
74 /* tridiagonal form. */
79 /* COMPZ (input) CHARACTER*1 */
80 /* = 'N': Compute eigenvalues only. */
81 /* = 'V': Compute eigenvalues and eigenvectors of the original */
82 /* symmetric matrix. On entry, Z must contain the */
83 /* orthogonal matrix used to reduce the original matrix */
84 /* to tridiagonal form. */
85 /* = 'I': Compute eigenvalues and eigenvectors of the */
86 /* tridiagonal matrix. Z is initialized to the identity */
89 /* N (input) INTEGER */
90 /* The order of the matrix. N >= 0. */
92 /* D (input/output) REAL array, dimension (N) */
93 /* On entry, the diagonal elements of the tridiagonal matrix. */
94 /* On exit, if INFO = 0, the eigenvalues in ascending order. */
96 /* E (input/output) REAL array, dimension (N-1) */
97 /* On entry, the (n-1) subdiagonal elements of the tridiagonal */
99 /* On exit, E has been destroyed. */
101 /* Z (input/output) REAL array, dimension (LDZ, N) */
102 /* On entry, if COMPZ = 'V', then Z contains the orthogonal */
103 /* matrix used in the reduction to tridiagonal form. */
104 /* On exit, if INFO = 0, then if COMPZ = 'V', Z contains the */
105 /* orthonormal eigenvectors of the original symmetric matrix, */
106 /* and if COMPZ = 'I', Z contains the orthonormal eigenvectors */
107 /* of the symmetric tridiagonal matrix. */
108 /* If COMPZ = 'N', then Z is not referenced. */
110 /* LDZ (input) INTEGER */
111 /* The leading dimension of the array Z. LDZ >= 1, and if */
112 /* eigenvectors are desired, then LDZ >= max(1,N). */
114 /* WORK (workspace) REAL array, dimension (max(1,2*N-2)) */
115 /* If COMPZ = 'N', then WORK is not referenced. */
117 /* INFO (output) INTEGER */
118 /* = 0: successful exit */
119 /* < 0: if INFO = -i, the i-th argument had an illegal value */
120 /* > 0: the algorithm has failed to find all the eigenvalues in */
121 /* a total of 30*N iterations; if INFO = i, then i */
122 /* elements of E have not converged to zero; on exit, D */
123 /* and E contain the elements of a symmetric tridiagonal */
124 /* matrix which is orthogonally similar to the original */
127 /* ===================================================================== */
129 /* .. Parameters .. */
131 /* .. Local Scalars .. */
133 /* .. External Functions .. */
135 /* .. External Subroutines .. */
137 /* .. Intrinsic Functions .. */
139 /* .. Executable Statements .. */
141 /* Test the input parameters. */
143 /* Parameter adjustments */
147 z_offset = 1 + z_dim1;
154 if (lsame_(compz, "N")) {
156 } else if (lsame_(compz, "V")) {
158 } else if (lsame_(compz, "I")) {
167 } else if (*ldz < 1 || icompz > 0 && *ldz < max(1,*n)) {
172 xerbla_("SSTEQR", &i__1);
176 /* Quick return if possible */
184 z__[z_dim1 + 1] = 1.f;
189 /* Determine the unit roundoff and over/underflow thresholds. */
192 /* Computing 2nd power */
195 safmin = slamch_("S");
196 safmax = 1.f / safmin;
197 ssfmax = sqrt(safmax) / 3.f;
198 ssfmin = sqrt(safmin) / eps2;
200 /* Compute the eigenvalues and eigenvectors of the tridiagonal */
204 slaset_("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
210 /* Determine where the matrix splits and choose QL or QR iteration */
211 /* for each block, according to whether top or bottom diagonal */
212 /* element is smaller. */
226 for (m = l1; m <= i__1; ++m) {
227 tst = (r__1 = e[m], dabs(r__1));
231 if (tst <= sqrt((r__1 = d__[m], dabs(r__1))) * sqrt((r__2 = d__[m
232 + 1], dabs(r__2))) * eps) {
251 /* Scale submatrix in rows and columns L to LEND */
254 anorm = slanst_("I", &i__1, &d__[l], &e[l]);
259 if (anorm > ssfmax) {
262 slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
265 slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
267 } else if (anorm < ssfmin) {
270 slascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
273 slascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
277 /* Choose between QL and QR iteration */
279 if ((r__1 = d__[lend], dabs(r__1)) < (r__2 = d__[l], dabs(r__2))) {
288 /* Look for small subdiagonal element. */
294 for (m = l; m <= i__1; ++m) {
295 /* Computing 2nd power */
296 r__2 = (r__1 = e[m], dabs(r__1));
298 if (tst <= eps2 * (r__1 = d__[m], dabs(r__1)) * (r__2 = d__[m
299 + 1], dabs(r__2)) + safmin) {
317 /* If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 */
318 /* to compute its eigensystem. */
322 slaev2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
324 work[*n - 1 + l] = s;
325 slasr_("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
326 z__[l * z_dim1 + 1], ldz);
328 slae2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
340 if (jtot == nmaxit) {
347 g = (d__[l + 1] - p) / (e[l] * 2.f);
348 r__ = slapy2_(&g, &c_b10);
349 g = d__[m] - p + e[l] / (g + r_sign(&r__, &g));
359 for (i__ = mm1; i__ >= i__1; --i__) {
362 slartg_(&g, &f, &c__, &s, &r__);
366 g = d__[i__ + 1] - p;
367 r__ = (d__[i__] - g) * s + c__ * 2.f * b;
369 d__[i__ + 1] = g + p;
372 /* If eigenvectors are desired, then save rotations. */
376 work[*n - 1 + i__] = -s;
382 /* If eigenvectors are desired, then apply saved rotations. */
386 slasr_("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l
394 /* Eigenvalue found. */
409 /* Look for small superdiagonal element. */
415 for (m = l; m >= i__1; --m) {
416 /* Computing 2nd power */
417 r__2 = (r__1 = e[m - 1], dabs(r__1));
419 if (tst <= eps2 * (r__1 = d__[m], dabs(r__1)) * (r__2 = d__[m
420 - 1], dabs(r__2)) + safmin) {
438 /* If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 */
439 /* to compute its eigensystem. */
443 slaev2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
446 work[*n - 1 + m] = s;
447 slasr_("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
448 z__[(l - 1) * z_dim1 + 1], ldz);
450 slae2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
462 if (jtot == nmaxit) {
469 g = (d__[l - 1] - p) / (e[l - 1] * 2.f);
470 r__ = slapy2_(&g, &c_b10);
471 g = d__[m] - p + e[l - 1] / (g + r_sign(&r__, &g));
481 for (i__ = m; i__ <= i__1; ++i__) {
484 slartg_(&g, &f, &c__, &s, &r__);
489 r__ = (d__[i__ + 1] - g) * s + c__ * 2.f * b;
494 /* If eigenvectors are desired, then save rotations. */
498 work[*n - 1 + i__] = s;
504 /* If eigenvectors are desired, then apply saved rotations. */
508 slasr_("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m
516 /* Eigenvalue found. */
529 /* Undo scaling if necessary */
533 i__1 = lendsv - lsv + 1;
534 slascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
537 slascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
539 } else if (iscale == 2) {
540 i__1 = lendsv - lsv + 1;
541 slascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
544 slascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
548 /* Check for no convergence to an eigenvalue after a total */
549 /* of N*MAXIT iterations. */
555 for (i__ = 1; i__ <= i__1; ++i__) {
563 /* Order eigenvalues and eigenvectors. */
570 slasrt_("I", n, &d__[1], info);
574 /* Use Selection Sort to minimize swaps of eigenvectors */
577 for (ii = 2; ii <= i__1; ++ii) {
582 for (j = ii; j <= i__2; ++j) {
592 sswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1],