3 /* Table of constant values */
5 static integer c__0 = 0;
6 static integer c__1 = 1;
7 static real c_b32 = 1.f;
9 /* Subroutine */ int ssterf_(integer *n, real *d__, real *e, integer *info)
11 /* System generated locals */
13 real r__1, r__2, r__3;
15 /* Builtin functions */
16 double sqrt(doublereal), r_sign(real *, real *);
23 real bb, rt1, rt2, eps, rte;
27 extern /* Subroutine */ int slae2_(real *, real *, real *, real *, real *)
29 real gamma, alpha, sigma, anorm;
30 extern doublereal slapy2_(real *, real *);
33 extern doublereal slamch_(char *);
35 extern /* Subroutine */ int xerbla_(char *, integer *);
37 extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
38 real *, integer *, integer *, real *, integer *, integer *);
43 extern doublereal slanst_(char *, integer *, real *, real *);
44 extern /* Subroutine */ int slasrt_(char *, integer *, real *, integer *);
47 /* -- LAPACK routine (version 3.1) -- */
48 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
51 /* .. Scalar Arguments .. */
53 /* .. Array Arguments .. */
59 /* SSTERF computes all eigenvalues of a symmetric tridiagonal matrix */
60 /* using the Pal-Walker-Kahan variant of the QL or QR algorithm. */
65 /* N (input) INTEGER */
66 /* The order of the matrix. N >= 0. */
68 /* D (input/output) REAL array, dimension (N) */
69 /* On entry, the n diagonal elements of the tridiagonal matrix. */
70 /* On exit, if INFO = 0, the eigenvalues in ascending order. */
72 /* E (input/output) REAL array, dimension (N-1) */
73 /* On entry, the (n-1) subdiagonal elements of the tridiagonal */
75 /* On exit, E has been destroyed. */
77 /* INFO (output) INTEGER */
78 /* = 0: successful exit */
79 /* < 0: if INFO = -i, the i-th argument had an illegal value */
80 /* > 0: the algorithm failed to find all of the eigenvalues in */
81 /* a total of 30*N iterations; if INFO = i, then i */
82 /* elements of E have not converged to zero. */
84 /* ===================================================================== */
86 /* .. Parameters .. */
88 /* .. Local Scalars .. */
90 /* .. External Functions .. */
92 /* .. External Subroutines .. */
94 /* .. Intrinsic Functions .. */
96 /* .. Executable Statements .. */
98 /* Test the input parameters. */
100 /* Parameter adjustments */
107 /* Quick return if possible */
112 xerbla_("SSTERF", &i__1);
119 /* Determine the unit roundoff for this environment. */
122 /* Computing 2nd power */
125 safmin = slamch_("S");
126 safmax = 1.f / safmin;
127 ssfmax = sqrt(safmax) / 3.f;
128 ssfmin = sqrt(safmin) / eps2;
130 /* Compute the eigenvalues of the tridiagonal matrix. */
136 /* Determine where the matrix splits and choose QL or QR iteration */
137 /* for each block, according to whether top or bottom diagonal */
138 /* element is smaller. */
150 for (m = l1; m <= i__1; ++m) {
151 if ((r__3 = e[m], dabs(r__3)) <= sqrt((r__1 = d__[m], dabs(r__1))) *
152 sqrt((r__2 = d__[m + 1], dabs(r__2))) * eps) {
170 /* Scale submatrix in rows and columns L to LEND */
173 anorm = slanst_("I", &i__1, &d__[l], &e[l]);
175 if (anorm > ssfmax) {
178 slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
181 slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
183 } else if (anorm < ssfmin) {
186 slascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
189 slascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
194 for (i__ = l; i__ <= i__1; ++i__) {
195 /* Computing 2nd power */
197 e[i__] = r__1 * r__1;
201 /* Choose between QL and QR iteration */
203 if ((r__1 = d__[lend], dabs(r__1)) < (r__2 = d__[l], dabs(r__2))) {
212 /* Look for small subdiagonal element. */
217 for (m = l; m <= i__1; ++m) {
218 if ((r__2 = e[m], dabs(r__2)) <= eps2 * (r__1 = d__[m] * d__[
219 m + 1], dabs(r__1))) {
236 /* If remaining matrix is 2 by 2, use SLAE2 to compute its */
241 slae2_(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
252 if (jtot == nmaxit) {
260 sigma = (d__[l + 1] - p) / (rte * 2.f);
261 r__ = slapy2_(&sigma, &c_b32);
262 sigma = p - rte / (sigma + r_sign(&r__, &sigma));
266 gamma = d__[m] - sigma;
272 for (i__ = m - 1; i__ >= i__1; --i__) {
276 e[i__ + 1] = s * r__;
283 gamma = c__ * (alpha - sigma) - s * oldgam;
284 d__[i__ + 1] = oldgam + (alpha - gamma);
286 p = gamma * gamma / c__;
294 d__[l] = sigma + gamma;
297 /* Eigenvalue found. */
312 /* Look for small superdiagonal element. */
316 for (m = l; m >= i__1; --m) {
317 if ((r__2 = e[m - 1], dabs(r__2)) <= eps2 * (r__1 = d__[m] * d__[
318 m - 1], dabs(r__1))) {
334 /* If remaining matrix is 2 by 2, use SLAE2 to compute its */
338 rte = sqrt(e[l - 1]);
339 slae2_(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
350 if (jtot == nmaxit) {
357 rte = sqrt(e[l - 1]);
358 sigma = (d__[l - 1] - p) / (rte * 2.f);
359 r__ = slapy2_(&sigma, &c_b32);
360 sigma = p - rte / (sigma + r_sign(&r__, &sigma));
364 gamma = d__[m] - sigma;
370 for (i__ = m; i__ <= i__1; ++i__) {
374 e[i__ - 1] = s * r__;
380 alpha = d__[i__ + 1];
381 gamma = c__ * (alpha - sigma) - s * oldgam;
382 d__[i__] = oldgam + (alpha - gamma);
384 p = gamma * gamma / c__;
392 d__[l] = sigma + gamma;
395 /* Eigenvalue found. */
408 /* Undo scaling if necessary */
412 i__1 = lendsv - lsv + 1;
413 slascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
417 i__1 = lendsv - lsv + 1;
418 slascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
422 /* Check for no convergence to an eigenvalue after a total */
423 /* of N*MAXIT iterations. */
429 for (i__ = 1; i__ <= i__1; ++i__) {
437 /* Sort eigenvalues in increasing order. */
440 slasrt_("I", n, &d__[1], info);