/* tql2.f -- translated by f2c (version 19961017).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Table of constant values */
static doublereal c_b10 = 1.;
/* Subroutine */ int tql2_(integer *nm, integer *n, doublereal *d__,
doublereal *e, doublereal *z__, integer *ierr)
{
/* System generated locals */
integer z_dim1, z_offset, i__1, i__2, i__3;
doublereal d__1, d__2;
/* Builtin functions */
double d_sign(doublereal *, doublereal *);
/* Local variables */
static doublereal c__, f, g, h__;
static integer i__, j, k, l, m;
static doublereal p, r__, s, c2, c3;
static integer l1, l2;
static doublereal s2;
static integer ii;
extern doublereal pythag_(doublereal *, doublereal *);
static doublereal dl1, el1;
static integer mml;
static doublereal tst1, tst2;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, */
/* NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND */
/* WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). */
/* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
/* OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. */
/* THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO */
/* BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS */
/* FULL MATRIX TO TRIDIAGONAL FORM. */
/* ON INPUT */
/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
/* DIMENSION STATEMENT. */
/* N IS THE ORDER OF THE MATRIX. */
/* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
/* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
/* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */
/* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */
/* REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS */
/* OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN */
/* THE IDENTITY MATRIX. */
/* ON OUTPUT */
/* D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN */
/* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT */
/* UNORDERED FOR INDICES 1,2,...,IERR-1. */
/* E HAS BEEN DESTROYED. */
/* Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC */
/* TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, */
/* Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED */
/* EIGENVALUES. */
/* IERR IS SET TO */
/* ZERO FOR NORMAL RETURN, */
/* J IF THE J-TH EIGENVALUE HAS NOT BEEN */
/* DETERMINED AFTER 30 ITERATIONS. */
/* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */
/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
*/
/* THIS VERSION DATED AUGUST 1983. */
/* ------------------------------------------------------------------
*/
/* Parameter adjustments */
z_dim1 = *nm;
z_offset = z_dim1 + 1;
z__ -= z_offset;
--e;
--d__;
/* Function Body */
*ierr = 0;
if (*n == 1) {
goto L1001;
}
i__1 = *n;
for (i__ = 2; i__ <= i__1; ++i__) {
/* L100: */
e[i__ - 1] = e[i__];
}
f = 0.;
tst1 = 0.;
e[*n] = 0.;
i__1 = *n;
for (l = 1; l <= i__1; ++l) {
j = 0;
h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2));
if (tst1 < h__) {
tst1 = h__;
}
/* .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
i__2 = *n;
for (m = l; m <= i__2; ++m) {
tst2 = tst1 + (d__1 = e[m], abs(d__1));
if (tst2 == tst1) {
goto L120;
}
/* .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT */
/* THROUGH THE BOTTOM OF THE LOOP .......... */
/* L110: */
}
L120:
if (m == l) {
goto L220;
}
L130:
if (j == 30) {
goto L1000;
}
++j;
/* .......... FORM SHIFT .......... */
l1 = l + 1;
l2 = l1 + 1;
g = d__[l];
p = (d__[l1] - g) / (e[l] * 2.);
r__ = pythag_(&p, &c_b10);
d__[l] = e[l] / (p + d_sign(&r__, &p));
d__[l1] = e[l] * (p + d_sign(&r__, &p));
dl1 = d__[l1];
h__ = g - d__[l];
if (l2 > *n) {
goto L145;
}
i__2 = *n;
for (i__ = l2; i__ <= i__2; ++i__) {
/* L140: */
d__[i__] -= h__;
}
L145:
f += h__;
/* .......... QL TRANSFORMATION .......... */
p = d__[m];
c__ = 1.;
c2 = c__;
el1 = e[l1];
s = 0.;
mml = m - l;
/* .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
i__2 = mml;
for (ii = 1; ii <= i__2; ++ii) {
c3 = c2;
c2 = c__;
s2 = s;
i__ = m - ii;
g = c__ * e[i__];
h__ = c__ * p;
r__ = pythag_(&p, &e[i__]);
e[i__ + 1] = s * r__;
s = e[i__] / r__;
c__ = p / r__;
p = c__ * d__[i__] - s * g;
d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]);
/* .......... FORM VECTOR .......... */
i__3 = *n;
for (k = 1; k <= i__3; ++k) {
h__ = z__[k + (i__ + 1) * z_dim1];
z__[k + (i__ + 1) * z_dim1] = s * z__[k + i__ * z_dim1] + c__
* h__;
z__[k + i__ * z_dim1] = c__ * z__[k + i__ * z_dim1] - s * h__;
/* L180: */
}
/* L200: */
}
p = -s * s2 * c3 * el1 * e[l] / dl1;
e[l] = s * p;
d__[l] = c__ * p;
tst2 = tst1 + (d__1 = e[l], abs(d__1));
if (tst2 > tst1) {
goto L130;
}
L220:
d__[l] += f;
/* L240: */
}
/* .......... ORDER EIGENVALUES AND EIGENVECTORS .......... */
i__1 = *n;
for (ii = 2; ii <= i__1; ++ii) {
i__ = ii - 1;
k = i__;
p = d__[i__];
i__2 = *n;
for (j = ii; j <= i__2; ++j) {
if (d__[j] >= p) {
goto L260;
}
k = j;
p = d__[j];
L260:
;
}
if (k == i__) {
goto L300;
}
d__[k] = d__[i__];
d__[i__] = p;
i__2 = *n;
for (j = 1; j <= i__2; ++j) {
p = z__[j + i__ * z_dim1];
z__[j + i__ * z_dim1] = z__[j + k * z_dim1];
z__[j + k * z_dim1] = p;
/* L280: */
}
L300:
;
}
goto L1001;
/* .......... SET ERROR -- NO CONVERGENCE TO AN */
/* EIGENVALUE AFTER 30 ITERATIONS .......... */
L1000:
*ierr = l;
L1001:
return 0;
} /* tql2_ */
syntax highlighted by Code2HTML, v. 0.9.1