/* hqr.f -- translated by f2c (version 19961017).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Subroutine */ int hqr_(integer *nm, integer *n, integer *low, integer *igh,
doublereal *h__, doublereal *wr, doublereal *wi, integer *ierr)
{
/* System generated locals */
integer h_dim1, h_offset, i__1, i__2, i__3;
doublereal d__1, d__2;
/* Builtin functions */
double sqrt(doublereal), d_sign(doublereal *, doublereal *);
/* Local variables */
static doublereal norm;
static integer i__, j, k, l, m;
static doublereal p, q, r__, s, t, w, x, y;
static integer na, en, ll, mm;
static doublereal zz;
static logical notlas;
static integer mp2, itn, its, enm2;
static doublereal tst1, tst2;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, */
/* NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). */
/* THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL */
/* UPPER HESSENBERG MATRIX BY THE QR METHOD. */
/* 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. */
/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
/* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
/* SET LOW=1, IGH=N. */
/* H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT */
/* THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG */
/* FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED */
/* IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */
/* ON OUTPUT */
/* H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED */
/* BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND */
/* BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. */
/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES */
/* ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */
/* OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */
/* HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN */
/* ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
/* FOR INDICES IERR+1,...,N. */
/* IERR IS SET TO */
/* ZERO FOR NORMAL RETURN, */
/* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
/* 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 */
--wi;
--wr;
h_dim1 = *nm;
h_offset = h_dim1 + 1;
h__ -= h_offset;
/* Function Body */
*ierr = 0;
norm = 0.;
k = 1;
/* .......... STORE ROOTS ISOLATED BY BALANC */
/* AND COMPUTE MATRIX NORM .......... */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = *n;
for (j = k; j <= i__2; ++j) {
/* L40: */
norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1));
}
k = i__;
if (i__ >= *low && i__ <= *igh) {
goto L50;
}
wr[i__] = h__[i__ + i__ * h_dim1];
wi[i__] = 0.;
L50:
;
}
en = *igh;
t = 0.;
itn = *n * 30;
/* .......... SEARCH FOR NEXT EIGENVALUES .......... */
L60:
if (en < *low) {
goto L1001;
}
its = 0;
na = en - 1;
enm2 = na - 1;
/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
/* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
L70:
i__1 = en;
for (ll = *low; ll <= i__1; ++ll) {
l = en + *low - ll;
if (l == *low) {
goto L100;
}
s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l
+ l * h_dim1], abs(d__2));
if (s == 0.) {
s = norm;
}
tst1 = s;
tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1));
if (tst2 == tst1) {
goto L100;
}
/* L80: */
}
/* .......... FORM SHIFT .......... */
L100:
x = h__[en + en * h_dim1];
if (l == en) {
goto L270;
}
y = h__[na + na * h_dim1];
w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
if (l == na) {
goto L280;
}
if (itn == 0) {
goto L1000;
}
if (its != 10 && its != 20) {
goto L130;
}
/* .......... FORM EXCEPTIONAL SHIFT .......... */
t += x;
i__1 = en;
for (i__ = *low; i__ <= i__1; ++i__) {
/* L120: */
h__[i__ + i__ * h_dim1] -= x;
}
s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 *
h_dim1], abs(d__2));
x = s * .75;
y = x;
w = s * -.4375 * s;
L130:
++its;
--itn;
/* .......... LOOK FOR TWO CONSECUTIVE SMALL */
/* SUB-DIAGONAL ELEMENTS. */
/* FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */
i__1 = enm2;
for (mm = l; mm <= i__1; ++mm) {
m = enm2 + l - mm;
zz = h__[m + m * h_dim1];
r__ = x - zz;
s = y - zz;
p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) *
h_dim1];
q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
r__ = h__[m + 2 + (m + 1) * h_dim1];
s = abs(p) + abs(q) + abs(r__);
p /= s;
q /= s;
r__ /= s;
if (m == l) {
goto L150;
}
tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) +
abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q)
+ abs(r__));
if (tst2 == tst1) {
goto L150;
}
/* L140: */
}
L150:
mp2 = m + 2;
i__1 = en;
for (i__ = mp2; i__ <= i__1; ++i__) {
h__[i__ + (i__ - 2) * h_dim1] = 0.;
if (i__ == mp2) {
goto L160;
}
h__[i__ + (i__ - 3) * h_dim1] = 0.;
L160:
;
}
/* .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */
/* COLUMNS M TO EN .......... */
i__1 = na;
for (k = m; k <= i__1; ++k) {
notlas = k != na;
if (k == m) {
goto L170;
}
p = h__[k + (k - 1) * h_dim1];
q = h__[k + 1 + (k - 1) * h_dim1];
r__ = 0.;
if (notlas) {
r__ = h__[k + 2 + (k - 1) * h_dim1];
}
x = abs(p) + abs(q) + abs(r__);
if (x == 0.) {
goto L260;
}
p /= x;
q /= x;
r__ /= x;
L170:
d__1 = sqrt(p * p + q * q + r__ * r__);
s = d_sign(&d__1, &p);
if (k == m) {
goto L180;
}
h__[k + (k - 1) * h_dim1] = -s * x;
goto L190;
L180:
if (l != m) {
h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
}
L190:
p += s;
x = p / s;
y = q / s;
zz = r__ / s;
q /= p;
r__ /= p;
if (notlas) {
goto L225;
}
/* .......... ROW MODIFICATION .......... */
i__2 = *n;
for (j = k; j <= i__2; ++j) {
p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
h__[k + j * h_dim1] -= p * x;
h__[k + 1 + j * h_dim1] -= p * y;
/* L200: */
}
/* Computing MIN */
i__2 = en, i__3 = k + 3;
j = min(i__2,i__3);
/* .......... COLUMN MODIFICATION .......... */
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
h__[i__ + k * h_dim1] -= p;
h__[i__ + (k + 1) * h_dim1] -= p * q;
/* L210: */
}
goto L255;
L225:
/* .......... ROW MODIFICATION .......... */
i__2 = *n;
for (j = k; j <= i__2; ++j) {
p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
k + 2 + j * h_dim1];
h__[k + j * h_dim1] -= p * x;
h__[k + 1 + j * h_dim1] -= p * y;
h__[k + 2 + j * h_dim1] -= p * zz;
/* L230: */
}
/* Computing MIN */
i__2 = en, i__3 = k + 3;
j = min(i__2,i__3);
/* .......... COLUMN MODIFICATION .......... */
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] +
zz * h__[i__ + (k + 2) * h_dim1];
h__[i__ + k * h_dim1] -= p;
h__[i__ + (k + 1) * h_dim1] -= p * q;
h__[i__ + (k + 2) * h_dim1] -= p * r__;
/* L240: */
}
L255:
L260:
;
}
goto L70;
/* .......... ONE ROOT FOUND .......... */
L270:
wr[en] = x + t;
wi[en] = 0.;
en = na;
goto L60;
/* .......... TWO ROOTS FOUND .......... */
L280:
p = (y - x) / 2.;
q = p * p + w;
zz = sqrt((abs(q)));
x += t;
if (q < 0.) {
goto L320;
}
/* .......... REAL PAIR .......... */
zz = p + d_sign(&zz, &p);
wr[na] = x + zz;
wr[en] = wr[na];
if (zz != 0.) {
wr[en] = x - w / zz;
}
wi[na] = 0.;
wi[en] = 0.;
goto L330;
/* .......... COMPLEX PAIR .......... */
L320:
wr[na] = x + p;
wr[en] = x + p;
wi[na] = zz;
wi[en] = -zz;
L330:
en = enm2;
goto L60;
/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
/* CONVERGED AFTER 30*N ITERATIONS .......... */
L1000:
*ierr = en;
L1001:
return 0;
} /* hqr_ */
syntax highlighted by Code2HTML, v. 0.9.1