/* comlr.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 comlr_(integer *nm, integer *n, integer *low, integer *
igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi,
integer *ierr)
{
/* System generated locals */
integer hr_dim1, hr_offset, hi_dim1, hi_offset, i__1, i__2;
doublereal d__1, d__2, d__3, d__4;
/* Local variables */
extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
, doublereal *, doublereal *, doublereal *);
static integer i__, j, l, m, en, ll, mm;
static doublereal si, ti, xi, yi, sr, tr, xr, yr;
static integer im1;
extern /* Subroutine */ int csroot_(doublereal *, doublereal *,
doublereal *, doublereal *);
static integer mp1, itn, its;
static doublereal zzi, zzr;
static integer enm1;
static doublereal tst1, tst2;
/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR, */
/* NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON. */
/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). */
/* THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX */
/* UPPER HESSENBERG MATRIX BY THE MODIFIED LR 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 CBAL. IF CBAL HAS NOT BEEN USED, */
/* SET LOW=1, IGH=N. */
/* HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
/* THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE */
/* MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY COMHES, */
/* IF PERFORMED. */
/* ON OUTPUT */
/* THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
/* DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE */
/* CALLING COMLR IF SUBSEQUENT CALCULATION OF */
/* EIGENVECTORS IS TO BE PERFORMED. */
/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
/* RESPECTIVELY, OF THE EIGENVALUES. 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. */
/* CALLS CDIV FOR COMPLEX DIVISION. */
/* CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
/* 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;
hi_dim1 = *nm;
hi_offset = hi_dim1 + 1;
hi -= hi_offset;
hr_dim1 = *nm;
hr_offset = hr_dim1 + 1;
hr -= hr_offset;
/* Function Body */
*ierr = 0;
/* .......... STORE ROOTS ISOLATED BY CBAL .......... */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
if (i__ >= *low && i__ <= *igh) {
goto L200;
}
wr[i__] = hr[i__ + i__ * hr_dim1];
wi[i__] = hi[i__ + i__ * hi_dim1];
L200:
;
}
en = *igh;
tr = 0.;
ti = 0.;
itn = *n * 30;
/* .......... SEARCH FOR NEXT EIGENVALUE .......... */
L220:
if (en < *low) {
goto L1001;
}
its = 0;
enm1 = en - 1;
/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
/* FOR L=EN STEP -1 UNTIL LOW D0 -- .......... */
L240:
i__1 = en;
for (ll = *low; ll <= i__1; ++ll) {
l = en + *low - ll;
if (l == *low) {
goto L300;
}
tst1 = (d__1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[
l - 1 + (l - 1) * hi_dim1], abs(d__2)) + (d__3 = hr[l + l *
hr_dim1], abs(d__3)) + (d__4 = hi[l + l * hi_dim1], abs(d__4))
;
tst2 = tst1 + (d__1 = hr[l + (l - 1) * hr_dim1], abs(d__1)) + (d__2 =
hi[l + (l - 1) * hi_dim1], abs(d__2));
if (tst2 == tst1) {
goto L300;
}
/* L260: */
}
/* .......... FORM SHIFT .......... */
L300:
if (l == en) {
goto L660;
}
if (itn == 0) {
goto L1000;
}
if (its == 10 || its == 20) {
goto L320;
}
sr = hr[en + en * hr_dim1];
si = hi[en + en * hi_dim1];
xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en *
hi_dim1] * hi[en + enm1 * hi_dim1];
xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en *
hi_dim1] * hr[en + enm1 * hr_dim1];
if (xr == 0. && xi == 0.) {
goto L340;
}
yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
/* Computing 2nd power */
d__2 = yr;
/* Computing 2nd power */
d__3 = yi;
d__1 = d__2 * d__2 - d__3 * d__3 + xr;
d__4 = yr * 2. * yi + xi;
csroot_(&d__1, &d__4, &zzr, &zzi);
if (yr * zzr + yi * zzi >= 0.) {
goto L310;
}
zzr = -zzr;
zzi = -zzi;
L310:
d__1 = yr + zzr;
d__2 = yi + zzi;
cdiv_(&xr, &xi, &d__1, &d__2, &xr, &xi);
sr -= xr;
si -= xi;
goto L340;
/* .......... FORM EXCEPTIONAL SHIFT .......... */
L320:
sr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hr[enm1 + (en
- 2) * hr_dim1], abs(d__2));
si = (d__1 = hi[en + enm1 * hi_dim1], abs(d__1)) + (d__2 = hi[enm1 + (en
- 2) * hi_dim1], abs(d__2));
L340:
i__1 = en;
for (i__ = *low; i__ <= i__1; ++i__) {
hr[i__ + i__ * hr_dim1] -= sr;
hi[i__ + i__ * hi_dim1] -= si;
/* L360: */
}
tr += sr;
ti += si;
++its;
--itn;
/* .......... LOOK FOR TWO CONSECUTIVE SMALL */
/* SUB-DIAGONAL ELEMENTS .......... */
xr = (d__1 = hr[enm1 + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[enm1 +
enm1 * hi_dim1], abs(d__2));
yr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[en + enm1 *
hi_dim1], abs(d__2));
zzr = (d__1 = hr[en + en * hr_dim1], abs(d__1)) + (d__2 = hi[en + en *
hi_dim1], abs(d__2));
/* .......... FOR M=EN-1 STEP -1 UNTIL L DO -- .......... */
i__1 = enm1;
for (mm = l; mm <= i__1; ++mm) {
m = enm1 + l - mm;
if (m == l) {
goto L420;
}
yi = yr;
yr = (d__1 = hr[m + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m + (
m - 1) * hi_dim1], abs(d__2));
xi = zzr;
zzr = xr;
xr = (d__1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m
- 1 + (m - 1) * hi_dim1], abs(d__2));
tst1 = zzr / yi * (zzr + xr + xi);
tst2 = tst1 + yr;
if (tst2 == tst1) {
goto L420;
}
/* L380: */
}
/* .......... TRIANGULAR DECOMPOSITION H=L*R .......... */
L420:
mp1 = m + 1;
i__1 = en;
for (i__ = mp1; i__ <= i__1; ++i__) {
im1 = i__ - 1;
xr = hr[im1 + im1 * hr_dim1];
xi = hi[im1 + im1 * hi_dim1];
yr = hr[i__ + im1 * hr_dim1];
yi = hi[i__ + im1 * hi_dim1];
if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
goto L460;
}
/* .......... INTERCHANGE ROWS OF HR AND HI .......... */
i__2 = en;
for (j = im1; j <= i__2; ++j) {
zzr = hr[im1 + j * hr_dim1];
hr[im1 + j * hr_dim1] = hr[i__ + j * hr_dim1];
hr[i__ + j * hr_dim1] = zzr;
zzi = hi[im1 + j * hi_dim1];
hi[im1 + j * hi_dim1] = hi[i__ + j * hi_dim1];
hi[i__ + j * hi_dim1] = zzi;
/* L440: */
}
cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
wr[i__] = 1.;
goto L480;
L460:
cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
wr[i__] = -1.;
L480:
hr[i__ + im1 * hr_dim1] = zzr;
hi[i__ + im1 * hi_dim1] = zzi;
i__2 = en;
for (j = i__; j <= i__2; ++j) {
hr[i__ + j * hr_dim1] = hr[i__ + j * hr_dim1] - zzr * hr[im1 + j *
hr_dim1] + zzi * hi[im1 + j * hi_dim1];
hi[i__ + j * hi_dim1] = hi[i__ + j * hi_dim1] - zzr * hi[im1 + j *
hi_dim1] - zzi * hr[im1 + j * hr_dim1];
/* L500: */
}
/* L520: */
}
/* .......... COMPOSITION R*L=H .......... */
i__1 = en;
for (j = mp1; j <= i__1; ++j) {
xr = hr[j + (j - 1) * hr_dim1];
xi = hi[j + (j - 1) * hi_dim1];
hr[j + (j - 1) * hr_dim1] = 0.;
hi[j + (j - 1) * hi_dim1] = 0.;
/* .......... INTERCHANGE COLUMNS OF HR AND HI, */
/* IF NECESSARY .......... */
if (wr[j] <= 0.) {
goto L580;
}
i__2 = j;
for (i__ = l; i__ <= i__2; ++i__) {
zzr = hr[i__ + (j - 1) * hr_dim1];
hr[i__ + (j - 1) * hr_dim1] = hr[i__ + j * hr_dim1];
hr[i__ + j * hr_dim1] = zzr;
zzi = hi[i__ + (j - 1) * hi_dim1];
hi[i__ + (j - 1) * hi_dim1] = hi[i__ + j * hi_dim1];
hi[i__ + j * hi_dim1] = zzi;
/* L540: */
}
L580:
i__2 = j;
for (i__ = l; i__ <= i__2; ++i__) {
hr[i__ + (j - 1) * hr_dim1] = hr[i__ + (j - 1) * hr_dim1] + xr *
hr[i__ + j * hr_dim1] - xi * hi[i__ + j * hi_dim1];
hi[i__ + (j - 1) * hi_dim1] = hi[i__ + (j - 1) * hi_dim1] + xr *
hi[i__ + j * hi_dim1] + xi * hr[i__ + j * hr_dim1];
/* L600: */
}
/* L640: */
}
goto L240;
/* .......... A ROOT FOUND .......... */
L660:
wr[en] = hr[en + en * hr_dim1] + tr;
wi[en] = hi[en + en * hi_dim1] + ti;
en = enm1;
goto L220;
/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
/* CONVERGED AFTER 30*N ITERATIONS .......... */
L1000:
*ierr = en;
L1001:
return 0;
} /* comlr_ */
syntax highlighted by Code2HTML, v. 0.9.1