/* imtql1.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 imtql1_(integer *n, doublereal *d__, doublereal *e,
integer *ierr)
{
/* System generated locals */
integer i__1, i__2;
doublereal d__1, d__2;
/* Builtin functions */
double d_sign(doublereal *, doublereal *);
/* Local variables */
static doublereal b, c__, f, g;
static integer i__, j, l, m;
static doublereal p, r__, s;
static integer ii;
extern doublereal pythag_(doublereal *, doublereal *);
static integer mml;
static doublereal tst1, tst2;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL1, */
/* NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, */
/* AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). */
/* THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC */
/* TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. */
/* ON INPUT */
/* 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. */
/* ON OUTPUT */
/* D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN */
/* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
/* ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
/* THE SMALLEST EIGENVALUES. */
/* E HAS BEEN DESTROYED. */
/* 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 */
--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__];
}
e[*n] = 0.;
i__1 = *n;
for (l = 1; l <= i__1; ++l) {
j = 0;
/* .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
L105:
i__2 = *n;
for (m = l; m <= i__2; ++m) {
if (m == *n) {
goto L120;
}
tst1 = (d__1 = d__[m], abs(d__1)) + (d__2 = d__[m + 1], abs(d__2))
;
tst2 = tst1 + (d__1 = e[m], abs(d__1));
if (tst2 == tst1) {
goto L120;
}
/* L110: */
}
L120:
p = d__[l];
if (m == l) {
goto L215;
}
if (j == 30) {
goto L1000;
}
++j;
/* .......... FORM SHIFT .......... */
g = (d__[l + 1] - p) / (e[l] * 2.);
r__ = pythag_(&g, &c_b10);
g = d__[m] - p + e[l] / (g + d_sign(&r__, &g));
s = 1.;
c__ = 1.;
p = 0.;
mml = m - l;
/* .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
i__2 = mml;
for (ii = 1; ii <= i__2; ++ii) {
i__ = m - ii;
f = s * e[i__];
b = c__ * e[i__];
r__ = pythag_(&f, &g);
e[i__ + 1] = r__;
if (r__ == 0.) {
goto L210;
}
s = f / r__;
c__ = g / r__;
g = d__[i__ + 1] - p;
r__ = (d__[i__] - g) * s + c__ * 2. * b;
p = s * r__;
d__[i__ + 1] = g + p;
g = c__ * r__ - b;
/* L200: */
}
d__[l] -= p;
e[l] = g;
e[m] = 0.;
goto L105;
/* .......... RECOVER FROM UNDERFLOW .......... */
L210:
d__[i__ + 1] -= p;
e[m] = 0.;
goto L105;
/* .......... ORDER EIGENVALUES .......... */
L215:
if (l == 1) {
goto L250;
}
/* .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
i__2 = l;
for (ii = 2; ii <= i__2; ++ii) {
i__ = l + 2 - ii;
if (p >= d__[i__ - 1]) {
goto L270;
}
d__[i__] = d__[i__ - 1];
/* L230: */
}
L250:
i__ = 1;
L270:
d__[i__] = p;
/* L290: */
}
goto L1001;
/* .......... SET ERROR -- NO CONVERGENCE TO AN */
/* EIGENVALUE AFTER 30 ITERATIONS .......... */
L1000:
*ierr = l;
L1001:
return 0;
} /* imtql1_ */
syntax highlighted by Code2HTML, v. 0.9.1