00001 /* hqr2.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 /* Table of constant values */
00009 
00010 static doublereal c_b49 = 0.;
00011 
00012 /* Subroutine */ int hqr2_(integer *nm, integer *n, integer *low, integer *
00013         igh, doublereal *h__, doublereal *wr, doublereal *wi, doublereal *z__,
00014          integer *ierr)
00015 {
00016     /* System generated locals */
00017     integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3;
00018     doublereal d__1, d__2, d__3, d__4;
00019 
00020     /* Builtin functions */
00021     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00022 
00023     /* Local variables */
00024     extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
00025             , doublereal *, doublereal *, doublereal *);
00026     static doublereal norm;
00027     static integer i__, j, k, l, m;
00028     static doublereal p, q, r__, s, t, w, x, y;
00029     static integer na, ii, en, jj;
00030     static doublereal ra, sa;
00031     static integer ll, mm, nn;
00032     static doublereal vi, vr, zz;
00033     static logical notlas;
00034     static integer mp2, itn, its, enm2;
00035     static doublereal tst1, tst2;
00036 
00037 
00038 
00039 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, */
00040 /*     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
00041 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
00042 
00043 /*     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
00044 /*     OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD.  THE */
00045 /*     EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND */
00046 /*     IF  ELMHES  AND  ELTRAN  OR  ORTHES  AND  ORTRAN  HAVE */
00047 /*     BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM */
00048 /*     AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS. */
00049 
00050 /*     ON INPUT */
00051 
00052 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00053 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00054 /*          DIMENSION STATEMENT. */
00055 
00056 /*        N IS THE ORDER OF THE MATRIX. */
00057 
00058 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00059 /*          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED, */
00060 /*          SET LOW=1, IGH=N. */
00061 
00062 /*        H CONTAINS THE UPPER HESSENBERG MATRIX. */
00063 
00064 /*        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY  ELTRAN */
00065 /*          AFTER THE REDUCTION BY  ELMHES, OR BY  ORTRAN  AFTER THE */
00066 /*          REDUCTION BY  ORTHES, IF PERFORMED.  IF THE EIGENVECTORS */
00067 /*          OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE */
00068 /*          IDENTITY MATRIX. */
00069 
00070 /*     ON OUTPUT */
00071 
00072 /*        H HAS BEEN DESTROYED. */
00073 
00074 /*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
00075 /*          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES */
00076 /*          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */
00077 /*          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */
00078 /*          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN */
00079 /*          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
00080 /*          FOR INDICES IERR+1,...,N. */
00081 
00082 /*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. */
00083 /*          IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z */
00084 /*          CONTAINS ITS EIGENVECTOR.  IF THE I-TH EIGENVALUE IS COMPLEX 
00085 */
00086 /*          WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH */
00087 /*          COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS */
00088 /*          EIGENVECTOR.  THE EIGENVECTORS ARE UNNORMALIZED.  IF AN */
00089 /*          ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND. 
00090 */
00091 
00092 /*        IERR IS SET TO */
00093 /*          ZERO       FOR NORMAL RETURN, */
00094 /*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
00095 /*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
00096 
00097 /*     CALLS CDIV FOR COMPLEX DIVISION. */
00098 
00099 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00100 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00101 */
00102 
00103 /*     THIS VERSION DATED AUGUST 1983. */
00104 
00105 /*     ------------------------------------------------------------------ 
00106 */
00107 
00108     /* Parameter adjustments */
00109     z_dim1 = *nm;
00110     z_offset = z_dim1 + 1;
00111     z__ -= z_offset;
00112     --wi;
00113     --wr;
00114     h_dim1 = *nm;
00115     h_offset = h_dim1 + 1;
00116     h__ -= h_offset;
00117 
00118     /* Function Body */
00119     *ierr = 0;
00120     norm = 0.;
00121     k = 1;
00122 /*     .......... STORE ROOTS ISOLATED BY BALANC */
00123 /*                AND COMPUTE MATRIX NORM .......... */
00124     i__1 = *n;
00125     for (i__ = 1; i__ <= i__1; ++i__) {
00126 
00127         i__2 = *n;
00128         for (j = k; j <= i__2; ++j) {
00129 /* L40: */
00130             norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1));
00131         }
00132 
00133         k = i__;
00134         if (i__ >= *low && i__ <= *igh) {
00135             goto L50;
00136         }
00137         wr[i__] = h__[i__ + i__ * h_dim1];
00138         wi[i__] = 0.;
00139 L50:
00140         ;
00141     }
00142 
00143     en = *igh;
00144     t = 0.;
00145     itn = *n * 30;
00146 /*     .......... SEARCH FOR NEXT EIGENVALUES .......... */
00147 L60:
00148     if (en < *low) {
00149         goto L340;
00150     }
00151     its = 0;
00152     na = en - 1;
00153     enm2 = na - 1;
00154 /*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
00155 /*                FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
00156 L70:
00157     i__1 = en;
00158     for (ll = *low; ll <= i__1; ++ll) {
00159         l = en + *low - ll;
00160         if (l == *low) {
00161             goto L100;
00162         }
00163         s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l 
00164                 + l * h_dim1], abs(d__2));
00165         if (s == 0.) {
00166             s = norm;
00167         }
00168         tst1 = s;
00169         tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1));
00170         if (tst2 == tst1) {
00171             goto L100;
00172         }
00173 /* L80: */
00174     }
00175 /*     .......... FORM SHIFT .......... */
00176 L100:
00177     x = h__[en + en * h_dim1];
00178     if (l == en) {
00179         goto L270;
00180     }
00181     y = h__[na + na * h_dim1];
00182     w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
00183     if (l == na) {
00184         goto L280;
00185     }
00186     if (itn == 0) {
00187         goto L1000;
00188     }
00189     if (its != 10 && its != 20) {
00190         goto L130;
00191     }
00192 /*     .......... FORM EXCEPTIONAL SHIFT .......... */
00193     t += x;
00194 
00195     i__1 = en;
00196     for (i__ = *low; i__ <= i__1; ++i__) {
00197 /* L120: */
00198         h__[i__ + i__ * h_dim1] -= x;
00199     }
00200 
00201     s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 * 
00202             h_dim1], abs(d__2));
00203     x = s * .75;
00204     y = x;
00205     w = s * -.4375 * s;
00206 L130:
00207     ++its;
00208     --itn;
00209 /*     .......... LOOK FOR TWO CONSECUTIVE SMALL */
00210 /*                SUB-DIAGONAL ELEMENTS. */
00211 /*                FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */
00212     i__1 = enm2;
00213     for (mm = l; mm <= i__1; ++mm) {
00214         m = enm2 + l - mm;
00215         zz = h__[m + m * h_dim1];
00216         r__ = x - zz;
00217         s = y - zz;
00218         p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) * 
00219                 h_dim1];
00220         q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
00221         r__ = h__[m + 2 + (m + 1) * h_dim1];
00222         s = abs(p) + abs(q) + abs(r__);
00223         p /= s;
00224         q /= s;
00225         r__ /= s;
00226         if (m == l) {
00227             goto L150;
00228         }
00229         tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) + 
00230                 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
00231         tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q) 
00232                 + abs(r__));
00233         if (tst2 == tst1) {
00234             goto L150;
00235         }
00236 /* L140: */
00237     }
00238 
00239 L150:
00240     mp2 = m + 2;
00241 
00242     i__1 = en;
00243     for (i__ = mp2; i__ <= i__1; ++i__) {
00244         h__[i__ + (i__ - 2) * h_dim1] = 0.;
00245         if (i__ == mp2) {
00246             goto L160;
00247         }
00248         h__[i__ + (i__ - 3) * h_dim1] = 0.;
00249 L160:
00250         ;
00251     }
00252 /*     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */
00253 /*                COLUMNS M TO EN .......... */
00254     i__1 = na;
00255     for (k = m; k <= i__1; ++k) {
00256         notlas = k != na;
00257         if (k == m) {
00258             goto L170;
00259         }
00260         p = h__[k + (k - 1) * h_dim1];
00261         q = h__[k + 1 + (k - 1) * h_dim1];
00262         r__ = 0.;
00263         if (notlas) {
00264             r__ = h__[k + 2 + (k - 1) * h_dim1];
00265         }
00266         x = abs(p) + abs(q) + abs(r__);
00267         if (x == 0.) {
00268             goto L260;
00269         }
00270         p /= x;
00271         q /= x;
00272         r__ /= x;
00273 L170:
00274         d__1 = sqrt(p * p + q * q + r__ * r__);
00275         s = d_sign(&d__1, &p);
00276         if (k == m) {
00277             goto L180;
00278         }
00279         h__[k + (k - 1) * h_dim1] = -s * x;
00280         goto L190;
00281 L180:
00282         if (l != m) {
00283             h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
00284         }
00285 L190:
00286         p += s;
00287         x = p / s;
00288         y = q / s;
00289         zz = r__ / s;
00290         q /= p;
00291         r__ /= p;
00292         if (notlas) {
00293             goto L225;
00294         }
00295 /*     .......... ROW MODIFICATION .......... */
00296         i__2 = *n;
00297         for (j = k; j <= i__2; ++j) {
00298             p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
00299             h__[k + j * h_dim1] -= p * x;
00300             h__[k + 1 + j * h_dim1] -= p * y;
00301 /* L200: */
00302         }
00303 
00304 /* Computing MIN */
00305         i__2 = en, i__3 = k + 3;
00306         j = min(i__2,i__3);
00307 /*     .......... COLUMN MODIFICATION .......... */
00308         i__2 = j;
00309         for (i__ = 1; i__ <= i__2; ++i__) {
00310             p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
00311             h__[i__ + k * h_dim1] -= p;
00312             h__[i__ + (k + 1) * h_dim1] -= p * q;
00313 /* L210: */
00314         }
00315 /*     .......... ACCUMULATE TRANSFORMATIONS .......... */
00316         i__2 = *igh;
00317         for (i__ = *low; i__ <= i__2; ++i__) {
00318             p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1];
00319             z__[i__ + k * z_dim1] -= p;
00320             z__[i__ + (k + 1) * z_dim1] -= p * q;
00321 /* L220: */
00322         }
00323         goto L255;
00324 L225:
00325 /*     .......... ROW MODIFICATION .......... */
00326         i__2 = *n;
00327         for (j = k; j <= i__2; ++j) {
00328             p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
00329                     k + 2 + j * h_dim1];
00330             h__[k + j * h_dim1] -= p * x;
00331             h__[k + 1 + j * h_dim1] -= p * y;
00332             h__[k + 2 + j * h_dim1] -= p * zz;
00333 /* L230: */
00334         }
00335 
00336 /* Computing MIN */
00337         i__2 = en, i__3 = k + 3;
00338         j = min(i__2,i__3);
00339 /*     .......... COLUMN MODIFICATION .......... */
00340         i__2 = j;
00341         for (i__ = 1; i__ <= i__2; ++i__) {
00342             p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] + 
00343                     zz * h__[i__ + (k + 2) * h_dim1];
00344             h__[i__ + k * h_dim1] -= p;
00345             h__[i__ + (k + 1) * h_dim1] -= p * q;
00346             h__[i__ + (k + 2) * h_dim1] -= p * r__;
00347 /* L240: */
00348         }
00349 /*     .......... ACCUMULATE TRANSFORMATIONS .......... */
00350         i__2 = *igh;
00351         for (i__ = *low; i__ <= i__2; ++i__) {
00352             p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1] + 
00353                     zz * z__[i__ + (k + 2) * z_dim1];
00354             z__[i__ + k * z_dim1] -= p;
00355             z__[i__ + (k + 1) * z_dim1] -= p * q;
00356             z__[i__ + (k + 2) * z_dim1] -= p * r__;
00357 /* L250: */
00358         }
00359 L255:
00360 
00361 L260:
00362         ;
00363     }
00364 
00365     goto L70;
00366 /*     .......... ONE ROOT FOUND .......... */
00367 L270:
00368     h__[en + en * h_dim1] = x + t;
00369     wr[en] = h__[en + en * h_dim1];
00370     wi[en] = 0.;
00371     en = na;
00372     goto L60;
00373 /*     .......... TWO ROOTS FOUND .......... */
00374 L280:
00375     p = (y - x) / 2.;
00376     q = p * p + w;
00377     zz = sqrt((abs(q)));
00378     h__[en + en * h_dim1] = x + t;
00379     x = h__[en + en * h_dim1];
00380     h__[na + na * h_dim1] = y + t;
00381     if (q < 0.) {
00382         goto L320;
00383     }
00384 /*     .......... REAL PAIR .......... */
00385     zz = p + d_sign(&zz, &p);
00386     wr[na] = x + zz;
00387     wr[en] = wr[na];
00388     if (zz != 0.) {
00389         wr[en] = x - w / zz;
00390     }
00391     wi[na] = 0.;
00392     wi[en] = 0.;
00393     x = h__[en + na * h_dim1];
00394     s = abs(x) + abs(zz);
00395     p = x / s;
00396     q = zz / s;
00397     r__ = sqrt(p * p + q * q);
00398     p /= r__;
00399     q /= r__;
00400 /*     .......... ROW MODIFICATION .......... */
00401     i__1 = *n;
00402     for (j = na; j <= i__1; ++j) {
00403         zz = h__[na + j * h_dim1];
00404         h__[na + j * h_dim1] = q * zz + p * h__[en + j * h_dim1];
00405         h__[en + j * h_dim1] = q * h__[en + j * h_dim1] - p * zz;
00406 /* L290: */
00407     }
00408 /*     .......... COLUMN MODIFICATION .......... */
00409     i__1 = en;
00410     for (i__ = 1; i__ <= i__1; ++i__) {
00411         zz = h__[i__ + na * h_dim1];
00412         h__[i__ + na * h_dim1] = q * zz + p * h__[i__ + en * h_dim1];
00413         h__[i__ + en * h_dim1] = q * h__[i__ + en * h_dim1] - p * zz;
00414 /* L300: */
00415     }
00416 /*     .......... ACCUMULATE TRANSFORMATIONS .......... */
00417     i__1 = *igh;
00418     for (i__ = *low; i__ <= i__1; ++i__) {
00419         zz = z__[i__ + na * z_dim1];
00420         z__[i__ + na * z_dim1] = q * zz + p * z__[i__ + en * z_dim1];
00421         z__[i__ + en * z_dim1] = q * z__[i__ + en * z_dim1] - p * zz;
00422 /* L310: */
00423     }
00424 
00425     goto L330;
00426 /*     .......... COMPLEX PAIR .......... */
00427 L320:
00428     wr[na] = x + p;
00429     wr[en] = x + p;
00430     wi[na] = zz;
00431     wi[en] = -zz;
00432 L330:
00433     en = enm2;
00434     goto L60;
00435 /*     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND */
00436 /*                VECTORS OF UPPER TRIANGULAR FORM .......... */
00437 L340:
00438     if (norm == 0.) {
00439         goto L1001;
00440     }
00441 /*     .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... */
00442     i__1 = *n;
00443     for (nn = 1; nn <= i__1; ++nn) {
00444         en = *n + 1 - nn;
00445         p = wr[en];
00446         q = wi[en];
00447         na = en - 1;
00448         if (q < 0.) {
00449             goto L710;
00450         } else if (q == 0) {
00451             goto L600;
00452         } else {
00453             goto L800;
00454         }
00455 /*     .......... REAL VECTOR .......... */
00456 L600:
00457         m = en;
00458         h__[en + en * h_dim1] = 1.;
00459         if (na == 0) {
00460             goto L800;
00461         }
00462 /*     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
00463         i__2 = na;
00464         for (ii = 1; ii <= i__2; ++ii) {
00465             i__ = en - ii;
00466             w = h__[i__ + i__ * h_dim1] - p;
00467             r__ = 0.;
00468 
00469             i__3 = en;
00470             for (j = m; j <= i__3; ++j) {
00471 /* L610: */
00472                 r__ += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
00473             }
00474 
00475             if (wi[i__] >= 0.) {
00476                 goto L630;
00477             }
00478             zz = w;
00479             s = r__;
00480             goto L700;
00481 L630:
00482             m = i__;
00483             if (wi[i__] != 0.) {
00484                 goto L640;
00485             }
00486             t = w;
00487             if (t != 0.) {
00488                 goto L635;
00489             }
00490             tst1 = norm;
00491             t = tst1;
00492 L632:
00493             t *= .01;
00494             tst2 = norm + t;
00495             if (tst2 > tst1) {
00496                 goto L632;
00497             }
00498 L635:
00499             h__[i__ + en * h_dim1] = -r__ / t;
00500             goto L680;
00501 /*     .......... SOLVE REAL EQUATIONS .......... */
00502 L640:
00503             x = h__[i__ + (i__ + 1) * h_dim1];
00504             y = h__[i__ + 1 + i__ * h_dim1];
00505             q = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__];
00506             t = (x * s - zz * r__) / q;
00507             h__[i__ + en * h_dim1] = t;
00508             if (abs(x) <= abs(zz)) {
00509                 goto L650;
00510             }
00511             h__[i__ + 1 + en * h_dim1] = (-r__ - w * t) / x;
00512             goto L680;
00513 L650:
00514             h__[i__ + 1 + en * h_dim1] = (-s - y * t) / zz;
00515 
00516 /*     .......... OVERFLOW CONTROL .......... */
00517 L680:
00518             t = (d__1 = h__[i__ + en * h_dim1], abs(d__1));
00519             if (t == 0.) {
00520                 goto L700;
00521             }
00522             tst1 = t;
00523             tst2 = tst1 + 1. / tst1;
00524             if (tst2 > tst1) {
00525                 goto L700;
00526             }
00527             i__3 = en;
00528             for (j = i__; j <= i__3; ++j) {
00529                 h__[j + en * h_dim1] /= t;
00530 /* L690: */
00531             }
00532 
00533 L700:
00534             ;
00535         }
00536 /*     .......... END REAL VECTOR .......... */
00537         goto L800;
00538 /*     .......... COMPLEX VECTOR .......... */
00539 L710:
00540         m = na;
00541 /*     .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT */
00542 /*                EIGENVECTOR MATRIX IS TRIANGULAR .......... */
00543         if ((d__1 = h__[en + na * h_dim1], abs(d__1)) <= (d__2 = h__[na + en *
00544                  h_dim1], abs(d__2))) {
00545             goto L720;
00546         }
00547         h__[na + na * h_dim1] = q / h__[en + na * h_dim1];
00548         h__[na + en * h_dim1] = -(h__[en + en * h_dim1] - p) / h__[en + na * 
00549                 h_dim1];
00550         goto L730;
00551 L720:
00552         d__1 = -h__[na + en * h_dim1];
00553         d__2 = h__[na + na * h_dim1] - p;
00554         cdiv_(&c_b49, &d__1, &d__2, &q, &h__[na + na * h_dim1], &h__[na + en *
00555                  h_dim1]);
00556 L730:
00557         h__[en + na * h_dim1] = 0.;
00558         h__[en + en * h_dim1] = 1.;
00559         enm2 = na - 1;
00560         if (enm2 == 0) {
00561             goto L800;
00562         }
00563 /*     .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... */
00564         i__2 = enm2;
00565         for (ii = 1; ii <= i__2; ++ii) {
00566             i__ = na - ii;
00567             w = h__[i__ + i__ * h_dim1] - p;
00568             ra = 0.;
00569             sa = 0.;
00570 
00571             i__3 = en;
00572             for (j = m; j <= i__3; ++j) {
00573                 ra += h__[i__ + j * h_dim1] * h__[j + na * h_dim1];
00574                 sa += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
00575 /* L760: */
00576             }
00577 
00578             if (wi[i__] >= 0.) {
00579                 goto L770;
00580             }
00581             zz = w;
00582             r__ = ra;
00583             s = sa;
00584             goto L795;
00585 L770:
00586             m = i__;
00587             if (wi[i__] != 0.) {
00588                 goto L780;
00589             }
00590             d__1 = -ra;
00591             d__2 = -sa;
00592             cdiv_(&d__1, &d__2, &w, &q, &h__[i__ + na * h_dim1], &h__[i__ + 
00593                     en * h_dim1]);
00594             goto L790;
00595 /*     .......... SOLVE COMPLEX EQUATIONS .......... */
00596 L780:
00597             x = h__[i__ + (i__ + 1) * h_dim1];
00598             y = h__[i__ + 1 + i__ * h_dim1];
00599             vr = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__] - q * q;
00600             vi = (wr[i__] - p) * 2. * q;
00601             if (vr != 0. || vi != 0.) {
00602                 goto L784;
00603             }
00604             tst1 = norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(zz));
00605             vr = tst1;
00606 L783:
00607             vr *= .01;
00608             tst2 = tst1 + vr;
00609             if (tst2 > tst1) {
00610                 goto L783;
00611             }
00612 L784:
00613             d__1 = x * r__ - zz * ra + q * sa;
00614             d__2 = x * s - zz * sa - q * ra;
00615             cdiv_(&d__1, &d__2, &vr, &vi, &h__[i__ + na * h_dim1], &h__[i__ + 
00616                     en * h_dim1]);
00617             if (abs(x) <= abs(zz) + abs(q)) {
00618                 goto L785;
00619             }
00620             h__[i__ + 1 + na * h_dim1] = (-ra - w * h__[i__ + na * h_dim1] + 
00621                     q * h__[i__ + en * h_dim1]) / x;
00622             h__[i__ + 1 + en * h_dim1] = (-sa - w * h__[i__ + en * h_dim1] - 
00623                     q * h__[i__ + na * h_dim1]) / x;
00624             goto L790;
00625 L785:
00626             d__1 = -r__ - y * h__[i__ + na * h_dim1];
00627             d__2 = -s - y * h__[i__ + en * h_dim1];
00628             cdiv_(&d__1, &d__2, &zz, &q, &h__[i__ + 1 + na * h_dim1], &h__[
00629                     i__ + 1 + en * h_dim1]);
00630 
00631 /*     .......... OVERFLOW CONTROL .......... */
00632 L790:
00633 /* Computing MAX */
00634             d__3 = (d__1 = h__[i__ + na * h_dim1], abs(d__1)), d__4 = (d__2 = 
00635                     h__[i__ + en * h_dim1], abs(d__2));
00636             t = max(d__3,d__4);
00637             if (t == 0.) {
00638                 goto L795;
00639             }
00640             tst1 = t;
00641             tst2 = tst1 + 1. / tst1;
00642             if (tst2 > tst1) {
00643                 goto L795;
00644             }
00645             i__3 = en;
00646             for (j = i__; j <= i__3; ++j) {
00647                 h__[j + na * h_dim1] /= t;
00648                 h__[j + en * h_dim1] /= t;
00649 /* L792: */
00650             }
00651 
00652 L795:
00653             ;
00654         }
00655 /*     .......... END COMPLEX VECTOR .......... */
00656 L800:
00657         ;
00658     }
00659 /*     .......... END BACK SUBSTITUTION. */
00660 /*                VECTORS OF ISOLATED ROOTS .......... */
00661     i__1 = *n;
00662     for (i__ = 1; i__ <= i__1; ++i__) {
00663         if (i__ >= *low && i__ <= *igh) {
00664             goto L840;
00665         }
00666 
00667         i__2 = *n;
00668         for (j = i__; j <= i__2; ++j) {
00669 /* L820: */
00670             z__[i__ + j * z_dim1] = h__[i__ + j * h_dim1];
00671         }
00672 
00673 L840:
00674         ;
00675     }
00676 /*     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
00677 /*                VECTORS OF ORIGINAL FULL MATRIX. */
00678 /*                FOR J=N STEP -1 UNTIL LOW DO -- .......... */
00679     i__1 = *n;
00680     for (jj = *low; jj <= i__1; ++jj) {
00681         j = *n + *low - jj;
00682         m = min(j,*igh);
00683 
00684         i__2 = *igh;
00685         for (i__ = *low; i__ <= i__2; ++i__) {
00686             zz = 0.;
00687 
00688             i__3 = m;
00689             for (k = *low; k <= i__3; ++k) {
00690 /* L860: */
00691                 zz += z__[i__ + k * z_dim1] * h__[k + j * h_dim1];
00692             }
00693 
00694             z__[i__ + j * z_dim1] = zz;
00695 /* L880: */
00696         }
00697     }
00698 
00699     goto L1001;
00700 /*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
00701 /*                CONVERGED AFTER 30*N ITERATIONS .......... */
00702 L1000:
00703     *ierr = en;
00704 L1001:
00705     return 0;
00706 } /* hqr2_ */
00707