00001 /* hqr.f -- translated by f2c (version 19961017).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 /* Subroutine */ int hqr_(integer *nm, integer *n, integer *low, integer *igh,
00009          doublereal *h__, doublereal *wr, doublereal *wi, integer *ierr)
00010 {
00011     /* System generated locals */
00012     integer h_dim1, h_offset, i__1, i__2, i__3;
00013     doublereal d__1, d__2;
00014 
00015     /* Builtin functions */
00016     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00017 
00018     /* Local variables */
00019     static doublereal norm;
00020     static integer i__, j, k, l, m;
00021     static doublereal p, q, r__, s, t, w, x, y;
00022     static integer na, en, ll, mm;
00023     static doublereal zz;
00024     static logical notlas;
00025     static integer mp2, itn, its, enm2;
00026     static doublereal tst1, tst2;
00027 
00028 
00029 
00030 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, */
00031 /*     NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. */
00032 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). */
00033 
00034 /*     THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL */
00035 /*     UPPER HESSENBERG MATRIX BY THE QR METHOD. */
00036 
00037 /*     ON INPUT */
00038 
00039 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00040 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00041 /*          DIMENSION STATEMENT. */
00042 
00043 /*        N IS THE ORDER OF THE MATRIX. */
00044 
00045 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00046 /*          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED, */
00047 /*          SET LOW=1, IGH=N. */
00048 
00049 /*        H CONTAINS THE UPPER HESSENBERG MATRIX.  INFORMATION ABOUT */
00050 /*          THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG */
00051 /*          FORM BY  ELMHES  OR  ORTHES, IF PERFORMED, IS STORED */
00052 /*          IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */
00053 
00054 /*     ON OUTPUT */
00055 
00056 /*        H HAS BEEN DESTROYED.  THEREFORE, IT MUST BE SAVED */
00057 /*          BEFORE CALLING  HQR  IF SUBSEQUENT CALCULATION AND */
00058 /*          BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. */
00059 
00060 /*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
00061 /*          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES */
00062 /*          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */
00063 /*          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */
00064 /*          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN */
00065 /*          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
00066 /*          FOR INDICES IERR+1,...,N. */
00067 
00068 /*        IERR IS SET TO */
00069 /*          ZERO       FOR NORMAL RETURN, */
00070 /*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
00071 /*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
00072 
00073 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00074 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00075 */
00076 
00077 /*     THIS VERSION DATED AUGUST 1983. */
00078 
00079 /*     ------------------------------------------------------------------ 
00080 */
00081 
00082     /* Parameter adjustments */
00083     --wi;
00084     --wr;
00085     h_dim1 = *nm;
00086     h_offset = h_dim1 + 1;
00087     h__ -= h_offset;
00088 
00089     /* Function Body */
00090     *ierr = 0;
00091     norm = 0.;
00092     k = 1;
00093 /*     .......... STORE ROOTS ISOLATED BY BALANC */
00094 /*                AND COMPUTE MATRIX NORM .......... */
00095     i__1 = *n;
00096     for (i__ = 1; i__ <= i__1; ++i__) {
00097 
00098         i__2 = *n;
00099         for (j = k; j <= i__2; ++j) {
00100 /* L40: */
00101             norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1));
00102         }
00103 
00104         k = i__;
00105         if (i__ >= *low && i__ <= *igh) {
00106             goto L50;
00107         }
00108         wr[i__] = h__[i__ + i__ * h_dim1];
00109         wi[i__] = 0.;
00110 L50:
00111         ;
00112     }
00113 
00114     en = *igh;
00115     t = 0.;
00116     itn = *n * 30;
00117 /*     .......... SEARCH FOR NEXT EIGENVALUES .......... */
00118 L60:
00119     if (en < *low) {
00120         goto L1001;
00121     }
00122     its = 0;
00123     na = en - 1;
00124     enm2 = na - 1;
00125 /*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
00126 /*                FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
00127 L70:
00128     i__1 = en;
00129     for (ll = *low; ll <= i__1; ++ll) {
00130         l = en + *low - ll;
00131         if (l == *low) {
00132             goto L100;
00133         }
00134         s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l 
00135                 + l * h_dim1], abs(d__2));
00136         if (s == 0.) {
00137             s = norm;
00138         }
00139         tst1 = s;
00140         tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1));
00141         if (tst2 == tst1) {
00142             goto L100;
00143         }
00144 /* L80: */
00145     }
00146 /*     .......... FORM SHIFT .......... */
00147 L100:
00148     x = h__[en + en * h_dim1];
00149     if (l == en) {
00150         goto L270;
00151     }
00152     y = h__[na + na * h_dim1];
00153     w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
00154     if (l == na) {
00155         goto L280;
00156     }
00157     if (itn == 0) {
00158         goto L1000;
00159     }
00160     if (its != 10 && its != 20) {
00161         goto L130;
00162     }
00163 /*     .......... FORM EXCEPTIONAL SHIFT .......... */
00164     t += x;
00165 
00166     i__1 = en;
00167     for (i__ = *low; i__ <= i__1; ++i__) {
00168 /* L120: */
00169         h__[i__ + i__ * h_dim1] -= x;
00170     }
00171 
00172     s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 * 
00173             h_dim1], abs(d__2));
00174     x = s * .75;
00175     y = x;
00176     w = s * -.4375 * s;
00177 L130:
00178     ++its;
00179     --itn;
00180 /*     .......... LOOK FOR TWO CONSECUTIVE SMALL */
00181 /*                SUB-DIAGONAL ELEMENTS. */
00182 /*                FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */
00183     i__1 = enm2;
00184     for (mm = l; mm <= i__1; ++mm) {
00185         m = enm2 + l - mm;
00186         zz = h__[m + m * h_dim1];
00187         r__ = x - zz;
00188         s = y - zz;
00189         p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) * 
00190                 h_dim1];
00191         q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
00192         r__ = h__[m + 2 + (m + 1) * h_dim1];
00193         s = abs(p) + abs(q) + abs(r__);
00194         p /= s;
00195         q /= s;
00196         r__ /= s;
00197         if (m == l) {
00198             goto L150;
00199         }
00200         tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) + 
00201                 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
00202         tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q) 
00203                 + abs(r__));
00204         if (tst2 == tst1) {
00205             goto L150;
00206         }
00207 /* L140: */
00208     }
00209 
00210 L150:
00211     mp2 = m + 2;
00212 
00213     i__1 = en;
00214     for (i__ = mp2; i__ <= i__1; ++i__) {
00215         h__[i__ + (i__ - 2) * h_dim1] = 0.;
00216         if (i__ == mp2) {
00217             goto L160;
00218         }
00219         h__[i__ + (i__ - 3) * h_dim1] = 0.;
00220 L160:
00221         ;
00222     }
00223 /*     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */
00224 /*                COLUMNS M TO EN .......... */
00225     i__1 = na;
00226     for (k = m; k <= i__1; ++k) {
00227         notlas = k != na;
00228         if (k == m) {
00229             goto L170;
00230         }
00231         p = h__[k + (k - 1) * h_dim1];
00232         q = h__[k + 1 + (k - 1) * h_dim1];
00233         r__ = 0.;
00234         if (notlas) {
00235             r__ = h__[k + 2 + (k - 1) * h_dim1];
00236         }
00237         x = abs(p) + abs(q) + abs(r__);
00238         if (x == 0.) {
00239             goto L260;
00240         }
00241         p /= x;
00242         q /= x;
00243         r__ /= x;
00244 L170:
00245         d__1 = sqrt(p * p + q * q + r__ * r__);
00246         s = d_sign(&d__1, &p);
00247         if (k == m) {
00248             goto L180;
00249         }
00250         h__[k + (k - 1) * h_dim1] = -s * x;
00251         goto L190;
00252 L180:
00253         if (l != m) {
00254             h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
00255         }
00256 L190:
00257         p += s;
00258         x = p / s;
00259         y = q / s;
00260         zz = r__ / s;
00261         q /= p;
00262         r__ /= p;
00263         if (notlas) {
00264             goto L225;
00265         }
00266 /*     .......... ROW MODIFICATION .......... */
00267         i__2 = *n;
00268         for (j = k; j <= i__2; ++j) {
00269             p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
00270             h__[k + j * h_dim1] -= p * x;
00271             h__[k + 1 + j * h_dim1] -= p * y;
00272 /* L200: */
00273         }
00274 
00275 /* Computing MIN */
00276         i__2 = en, i__3 = k + 3;
00277         j = min(i__2,i__3);
00278 /*     .......... COLUMN MODIFICATION .......... */
00279         i__2 = j;
00280         for (i__ = 1; i__ <= i__2; ++i__) {
00281             p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
00282             h__[i__ + k * h_dim1] -= p;
00283             h__[i__ + (k + 1) * h_dim1] -= p * q;
00284 /* L210: */
00285         }
00286         goto L255;
00287 L225:
00288 /*     .......... ROW MODIFICATION .......... */
00289         i__2 = *n;
00290         for (j = k; j <= i__2; ++j) {
00291             p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
00292                     k + 2 + j * h_dim1];
00293             h__[k + j * h_dim1] -= p * x;
00294             h__[k + 1 + j * h_dim1] -= p * y;
00295             h__[k + 2 + j * h_dim1] -= p * zz;
00296 /* L230: */
00297         }
00298 
00299 /* Computing MIN */
00300         i__2 = en, i__3 = k + 3;
00301         j = min(i__2,i__3);
00302 /*     .......... COLUMN MODIFICATION .......... */
00303         i__2 = j;
00304         for (i__ = 1; i__ <= i__2; ++i__) {
00305             p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] + 
00306                     zz * h__[i__ + (k + 2) * h_dim1];
00307             h__[i__ + k * h_dim1] -= p;
00308             h__[i__ + (k + 1) * h_dim1] -= p * q;
00309             h__[i__ + (k + 2) * h_dim1] -= p * r__;
00310 /* L240: */
00311         }
00312 L255:
00313 
00314 L260:
00315         ;
00316     }
00317 
00318     goto L70;
00319 /*     .......... ONE ROOT FOUND .......... */
00320 L270:
00321     wr[en] = x + t;
00322     wi[en] = 0.;
00323     en = na;
00324     goto L60;
00325 /*     .......... TWO ROOTS FOUND .......... */
00326 L280:
00327     p = (y - x) / 2.;
00328     q = p * p + w;
00329     zz = sqrt((abs(q)));
00330     x += t;
00331     if (q < 0.) {
00332         goto L320;
00333     }
00334 /*     .......... REAL PAIR .......... */
00335     zz = p + d_sign(&zz, &p);
00336     wr[na] = x + zz;
00337     wr[en] = wr[na];
00338     if (zz != 0.) {
00339         wr[en] = x - w / zz;
00340     }
00341     wi[na] = 0.;
00342     wi[en] = 0.;
00343     goto L330;
00344 /*     .......... COMPLEX PAIR .......... */
00345 L320:
00346     wr[na] = x + p;
00347     wr[en] = x + p;
00348     wi[na] = zz;
00349     wi[en] = -zz;
00350 L330:
00351     en = enm2;
00352     goto L60;
00353 /*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
00354 /*                CONVERGED AFTER 30*N ITERATIONS .......... */
00355 L1000:
00356     *ierr = en;
00357 L1001:
00358     return 0;
00359 } /* hqr_ */
00360