/* 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_ */