00001 /* comqr2.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 comqr2_(integer *nm, integer *n, integer *low, integer *
00009         igh, doublereal *ortr, doublereal *orti, doublereal *hr, doublereal *
00010         hi, doublereal *wr, doublereal *wi, doublereal *zr, doublereal *zi, 
00011         integer *ierr)
00012 {
00013     /* System generated locals */
00014     integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset, 
00015             zi_dim1, zi_offset, i__1, i__2, i__3;
00016     doublereal d__1, d__2, d__3, d__4;
00017 
00018     /* Local variables */
00019     static integer iend;
00020     extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
00021             , doublereal *, doublereal *, doublereal *);
00022     static doublereal norm;
00023     static integer i__, j, k, l, m, ii, en, jj, ll, nn;
00024     static doublereal si, ti, xi, yi, sr, tr, xr, yr;
00025     extern doublereal pythag_(doublereal *, doublereal *);
00026     extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 
00027             doublereal *, doublereal *);
00028     static integer ip1, lp1, itn, its;
00029     static doublereal zzi, zzr;
00030     static integer enm1;
00031     static doublereal tst1, tst2;
00032 
00033 
00034 
00035 /*     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE */
00036 /*     ALGOL PROCEDURE  COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS */
00037 /*     AND WILKINSON. */
00038 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
00039 /*     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS */
00040 /*     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. */
00041 
00042 /*     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
00043 /*     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR */
00044 /*     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX */
00045 /*     CAN ALSO BE FOUND IF  CORTH  HAS BEEN USED TO REDUCE */
00046 /*     THIS GENERAL MATRIX TO HESSENBERG FORM. */
00047 
00048 /*     ON INPUT */
00049 
00050 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00051 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00052 /*          DIMENSION STATEMENT. */
00053 
00054 /*        N IS THE ORDER OF THE MATRIX. */
00055 
00056 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00057 /*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
00058 /*          SET LOW=1, IGH=N. */
00059 
00060 /*        ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */
00061 /*          FORMATIONS USED IN THE REDUCTION BY  CORTH, IF PERFORMED. */
00062 /*          ONLY ELEMENTS LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS 
00063 */
00064 /*          OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND */
00065 /*          ORTI(J) TO 0.0D0 FOR THESE ELEMENTS. */
00066 
00067 /*        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
00068 /*          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
00069 /*          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER */
00070 /*          INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE 
00071 */
00072 /*          REDUCTION BY  CORTH, IF PERFORMED.  IF THE EIGENVECTORS OF */
00073 /*          THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE */
00074 /*          ARBITRARY. */
00075 
00076 /*     ON OUTPUT */
00077 
00078 /*        ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI */
00079 /*          HAVE BEEN DESTROYED. */
00080 
00081 /*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
00082 /*          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR */
00083 /*          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
00084 /*          FOR INDICES IERR+1,...,N. */
00085 
00086 /*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
00087 /*          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS */
00088 /*          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF */
00089 /*          THE EIGENVECTORS HAS BEEN FOUND. */
00090 
00091 /*        IERR IS SET TO */
00092 /*          ZERO       FOR NORMAL RETURN, */
00093 /*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
00094 /*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
00095 
00096 /*     CALLS CDIV FOR COMPLEX DIVISION. */
00097 /*     CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
00098 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00099 
00100 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00101 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00102 */
00103 
00104 /*     THIS VERSION DATED AUGUST 1983. */
00105 
00106 /*     ------------------------------------------------------------------ 
00107 */
00108 
00109     /* Parameter adjustments */
00110     zi_dim1 = *nm;
00111     zi_offset = zi_dim1 + 1;
00112     zi -= zi_offset;
00113     zr_dim1 = *nm;
00114     zr_offset = zr_dim1 + 1;
00115     zr -= zr_offset;
00116     --wi;
00117     --wr;
00118     hi_dim1 = *nm;
00119     hi_offset = hi_dim1 + 1;
00120     hi -= hi_offset;
00121     hr_dim1 = *nm;
00122     hr_offset = hr_dim1 + 1;
00123     hr -= hr_offset;
00124     --orti;
00125     --ortr;
00126 
00127     /* Function Body */
00128     *ierr = 0;
00129 /*     .......... INITIALIZE EIGENVECTOR MATRIX .......... */
00130     i__1 = *n;
00131     for (j = 1; j <= i__1; ++j) {
00132 
00133         i__2 = *n;
00134         for (i__ = 1; i__ <= i__2; ++i__) {
00135             zr[i__ + j * zr_dim1] = 0.;
00136             zi[i__ + j * zi_dim1] = 0.;
00137 /* L100: */
00138         }
00139         zr[j + j * zr_dim1] = 1.;
00140 /* L101: */
00141     }
00142 /*     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS */
00143 /*                FROM THE INFORMATION LEFT BY CORTH .......... */
00144     iend = *igh - *low - 1;
00145     if (iend < 0) {
00146         goto L180;
00147     } else if (iend == 0) {
00148         goto L150;
00149     } else {
00150         goto L105;
00151     }
00152 /*     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
00153 L105:
00154     i__1 = iend;
00155     for (ii = 1; ii <= i__1; ++ii) {
00156         i__ = *igh - ii;
00157         if (ortr[i__] == 0. && orti[i__] == 0.) {
00158             goto L140;
00159         }
00160         if (hr[i__ + (i__ - 1) * hr_dim1] == 0. && hi[i__ + (i__ - 1) * 
00161                 hi_dim1] == 0.) {
00162             goto L140;
00163         }
00164 /*     .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ........
00165 .. */
00166         norm = hr[i__ + (i__ - 1) * hr_dim1] * ortr[i__] + hi[i__ + (i__ - 1) 
00167                 * hi_dim1] * orti[i__];
00168         ip1 = i__ + 1;
00169 
00170         i__2 = *igh;
00171         for (k = ip1; k <= i__2; ++k) {
00172             ortr[k] = hr[k + (i__ - 1) * hr_dim1];
00173             orti[k] = hi[k + (i__ - 1) * hi_dim1];
00174 /* L110: */
00175         }
00176 
00177         i__2 = *igh;
00178         for (j = i__; j <= i__2; ++j) {
00179             sr = 0.;
00180             si = 0.;
00181 
00182             i__3 = *igh;
00183             for (k = i__; k <= i__3; ++k) {
00184                 sr = sr + ortr[k] * zr[k + j * zr_dim1] + orti[k] * zi[k + j *
00185                          zi_dim1];
00186                 si = si + ortr[k] * zi[k + j * zi_dim1] - orti[k] * zr[k + j *
00187                          zr_dim1];
00188 /* L115: */
00189             }
00190 
00191             sr /= norm;
00192             si /= norm;
00193 
00194             i__3 = *igh;
00195             for (k = i__; k <= i__3; ++k) {
00196                 zr[k + j * zr_dim1] = zr[k + j * zr_dim1] + sr * ortr[k] - si 
00197                         * orti[k];
00198                 zi[k + j * zi_dim1] = zi[k + j * zi_dim1] + sr * orti[k] + si 
00199                         * ortr[k];
00200 /* L120: */
00201             }
00202 
00203 /* L130: */
00204         }
00205 
00206 L140:
00207         ;
00208     }
00209 /*     .......... CREATE REAL SUBDIAGONAL ELEMENTS .......... */
00210 L150:
00211     l = *low + 1;
00212 
00213     i__1 = *igh;
00214     for (i__ = l; i__ <= i__1; ++i__) {
00215 /* Computing MIN */
00216         i__2 = i__ + 1;
00217         ll = min(i__2,*igh);
00218         if (hi[i__ + (i__ - 1) * hi_dim1] == 0.) {
00219             goto L170;
00220         }
00221         norm = pythag_(&hr[i__ + (i__ - 1) * hr_dim1], &hi[i__ + (i__ - 1) * 
00222                 hi_dim1]);
00223         yr = hr[i__ + (i__ - 1) * hr_dim1] / norm;
00224         yi = hi[i__ + (i__ - 1) * hi_dim1] / norm;
00225         hr[i__ + (i__ - 1) * hr_dim1] = norm;
00226         hi[i__ + (i__ - 1) * hi_dim1] = 0.;
00227 
00228         i__2 = *n;
00229         for (j = i__; j <= i__2; ++j) {
00230             si = yr * hi[i__ + j * hi_dim1] - yi * hr[i__ + j * hr_dim1];
00231             hr[i__ + j * hr_dim1] = yr * hr[i__ + j * hr_dim1] + yi * hi[i__ 
00232                     + j * hi_dim1];
00233             hi[i__ + j * hi_dim1] = si;
00234 /* L155: */
00235         }
00236 
00237         i__2 = ll;
00238         for (j = 1; j <= i__2; ++j) {
00239             si = yr * hi[j + i__ * hi_dim1] + yi * hr[j + i__ * hr_dim1];
00240             hr[j + i__ * hr_dim1] = yr * hr[j + i__ * hr_dim1] - yi * hi[j + 
00241                     i__ * hi_dim1];
00242             hi[j + i__ * hi_dim1] = si;
00243 /* L160: */
00244         }
00245 
00246         i__2 = *igh;
00247         for (j = *low; j <= i__2; ++j) {
00248             si = yr * zi[j + i__ * zi_dim1] + yi * zr[j + i__ * zr_dim1];
00249             zr[j + i__ * zr_dim1] = yr * zr[j + i__ * zr_dim1] - yi * zi[j + 
00250                     i__ * zi_dim1];
00251             zi[j + i__ * zi_dim1] = si;
00252 /* L165: */
00253         }
00254 
00255 L170:
00256         ;
00257     }
00258 /*     .......... STORE ROOTS ISOLATED BY CBAL .......... */
00259 L180:
00260     i__1 = *n;
00261     for (i__ = 1; i__ <= i__1; ++i__) {
00262         if (i__ >= *low && i__ <= *igh) {
00263             goto L200;
00264         }
00265         wr[i__] = hr[i__ + i__ * hr_dim1];
00266         wi[i__] = hi[i__ + i__ * hi_dim1];
00267 L200:
00268         ;
00269     }
00270 
00271     en = *igh;
00272     tr = 0.;
00273     ti = 0.;
00274     itn = *n * 30;
00275 /*     .......... SEARCH FOR NEXT EIGENVALUE .......... */
00276 L220:
00277     if (en < *low) {
00278         goto L680;
00279     }
00280     its = 0;
00281     enm1 = en - 1;
00282 /*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
00283 /*                FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
00284 L240:
00285     i__1 = en;
00286     for (ll = *low; ll <= i__1; ++ll) {
00287         l = en + *low - ll;
00288         if (l == *low) {
00289             goto L300;
00290         }
00291         tst1 = (d__1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[
00292                 l - 1 + (l - 1) * hi_dim1], abs(d__2)) + (d__3 = hr[l + l * 
00293                 hr_dim1], abs(d__3)) + (d__4 = hi[l + l * hi_dim1], abs(d__4))
00294                 ;
00295         tst2 = tst1 + (d__1 = hr[l + (l - 1) * hr_dim1], abs(d__1));
00296         if (tst2 == tst1) {
00297             goto L300;
00298         }
00299 /* L260: */
00300     }
00301 /*     .......... FORM SHIFT .......... */
00302 L300:
00303     if (l == en) {
00304         goto L660;
00305     }
00306     if (itn == 0) {
00307         goto L1000;
00308     }
00309     if (its == 10 || its == 20) {
00310         goto L320;
00311     }
00312     sr = hr[en + en * hr_dim1];
00313     si = hi[en + en * hi_dim1];
00314     xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
00315     xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
00316     if (xr == 0. && xi == 0.) {
00317         goto L340;
00318     }
00319     yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
00320     yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
00321 /* Computing 2nd power */
00322     d__2 = yr;
00323 /* Computing 2nd power */
00324     d__3 = yi;
00325     d__1 = d__2 * d__2 - d__3 * d__3 + xr;
00326     d__4 = yr * 2. * yi + xi;
00327     csroot_(&d__1, &d__4, &zzr, &zzi);
00328     if (yr * zzr + yi * zzi >= 0.) {
00329         goto L310;
00330     }
00331     zzr = -zzr;
00332     zzi = -zzi;
00333 L310:
00334     d__1 = yr + zzr;
00335     d__2 = yi + zzi;
00336     cdiv_(&xr, &xi, &d__1, &d__2, &xr, &xi);
00337     sr -= xr;
00338     si -= xi;
00339     goto L340;
00340 /*     .......... FORM EXCEPTIONAL SHIFT .......... */
00341 L320:
00342     sr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hr[enm1 + (en 
00343             - 2) * hr_dim1], abs(d__2));
00344     si = 0.;
00345 
00346 L340:
00347     i__1 = en;
00348     for (i__ = *low; i__ <= i__1; ++i__) {
00349         hr[i__ + i__ * hr_dim1] -= sr;
00350         hi[i__ + i__ * hi_dim1] -= si;
00351 /* L360: */
00352     }
00353 
00354     tr += sr;
00355     ti += si;
00356     ++its;
00357     --itn;
00358 /*     .......... REDUCE TO TRIANGLE (ROWS) .......... */
00359     lp1 = l + 1;
00360 
00361     i__1 = en;
00362     for (i__ = lp1; i__ <= i__1; ++i__) {
00363         sr = hr[i__ + (i__ - 1) * hr_dim1];
00364         hr[i__ + (i__ - 1) * hr_dim1] = 0.;
00365         d__1 = pythag_(&hr[i__ - 1 + (i__ - 1) * hr_dim1], &hi[i__ - 1 + (i__ 
00366                 - 1) * hi_dim1]);
00367         norm = pythag_(&d__1, &sr);
00368         xr = hr[i__ - 1 + (i__ - 1) * hr_dim1] / norm;
00369         wr[i__ - 1] = xr;
00370         xi = hi[i__ - 1 + (i__ - 1) * hi_dim1] / norm;
00371         wi[i__ - 1] = xi;
00372         hr[i__ - 1 + (i__ - 1) * hr_dim1] = norm;
00373         hi[i__ - 1 + (i__ - 1) * hi_dim1] = 0.;
00374         hi[i__ + (i__ - 1) * hi_dim1] = sr / norm;
00375 
00376         i__2 = *n;
00377         for (j = i__; j <= i__2; ++j) {
00378             yr = hr[i__ - 1 + j * hr_dim1];
00379             yi = hi[i__ - 1 + j * hi_dim1];
00380             zzr = hr[i__ + j * hr_dim1];
00381             zzi = hi[i__ + j * hi_dim1];
00382             hr[i__ - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i__ + (i__ - 1)
00383                      * hi_dim1] * zzr;
00384             hi[i__ - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i__ + (i__ - 1)
00385                      * hi_dim1] * zzi;
00386             hr[i__ + j * hr_dim1] = xr * zzr - xi * zzi - hi[i__ + (i__ - 1) *
00387                      hi_dim1] * yr;
00388             hi[i__ + j * hi_dim1] = xr * zzi + xi * zzr - hi[i__ + (i__ - 1) *
00389                      hi_dim1] * yi;
00390 /* L490: */
00391         }
00392 
00393 /* L500: */
00394     }
00395 
00396     si = hi[en + en * hi_dim1];
00397     if (si == 0.) {
00398         goto L540;
00399     }
00400     norm = pythag_(&hr[en + en * hr_dim1], &si);
00401     sr = hr[en + en * hr_dim1] / norm;
00402     si /= norm;
00403     hr[en + en * hr_dim1] = norm;
00404     hi[en + en * hi_dim1] = 0.;
00405     if (en == *n) {
00406         goto L540;
00407     }
00408     ip1 = en + 1;
00409 
00410     i__1 = *n;
00411     for (j = ip1; j <= i__1; ++j) {
00412         yr = hr[en + j * hr_dim1];
00413         yi = hi[en + j * hi_dim1];
00414         hr[en + j * hr_dim1] = sr * yr + si * yi;
00415         hi[en + j * hi_dim1] = sr * yi - si * yr;
00416 /* L520: */
00417     }
00418 /*     .......... INVERSE OPERATION (COLUMNS) .......... */
00419 L540:
00420     i__1 = en;
00421     for (j = lp1; j <= i__1; ++j) {
00422         xr = wr[j - 1];
00423         xi = wi[j - 1];
00424 
00425         i__2 = j;
00426         for (i__ = 1; i__ <= i__2; ++i__) {
00427             yr = hr[i__ + (j - 1) * hr_dim1];
00428             yi = 0.;
00429             zzr = hr[i__ + j * hr_dim1];
00430             zzi = hi[i__ + j * hi_dim1];
00431             if (i__ == j) {
00432                 goto L560;
00433             }
00434             yi = hi[i__ + (j - 1) * hi_dim1];
00435             hi[i__ + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
00436                      hi_dim1] * zzi;
00437 L560:
00438             hr[i__ + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
00439                      hi_dim1] * zzr;
00440             hr[i__ + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
00441                     hi_dim1] * yr;
00442             hi[i__ + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
00443                     hi_dim1] * yi;
00444 /* L580: */
00445         }
00446 
00447         i__2 = *igh;
00448         for (i__ = *low; i__ <= i__2; ++i__) {
00449             yr = zr[i__ + (j - 1) * zr_dim1];
00450             yi = zi[i__ + (j - 1) * zi_dim1];
00451             zzr = zr[i__ + j * zr_dim1];
00452             zzi = zi[i__ + j * zi_dim1];
00453             zr[i__ + (j - 1) * zr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
00454                      hi_dim1] * zzr;
00455             zi[i__ + (j - 1) * zi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
00456                      hi_dim1] * zzi;
00457             zr[i__ + j * zr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 
00458                     hi_dim1] * yr;
00459             zi[i__ + j * zi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 
00460                     hi_dim1] * yi;
00461 /* L590: */
00462         }
00463 
00464 /* L600: */
00465     }
00466 
00467     if (si == 0.) {
00468         goto L240;
00469     }
00470 
00471     i__1 = en;
00472     for (i__ = 1; i__ <= i__1; ++i__) {
00473         yr = hr[i__ + en * hr_dim1];
00474         yi = hi[i__ + en * hi_dim1];
00475         hr[i__ + en * hr_dim1] = sr * yr - si * yi;
00476         hi[i__ + en * hi_dim1] = sr * yi + si * yr;
00477 /* L630: */
00478     }
00479 
00480     i__1 = *igh;
00481     for (i__ = *low; i__ <= i__1; ++i__) {
00482         yr = zr[i__ + en * zr_dim1];
00483         yi = zi[i__ + en * zi_dim1];
00484         zr[i__ + en * zr_dim1] = sr * yr - si * yi;
00485         zi[i__ + en * zi_dim1] = sr * yi + si * yr;
00486 /* L640: */
00487     }
00488 
00489     goto L240;
00490 /*     .......... A ROOT FOUND .......... */
00491 L660:
00492     hr[en + en * hr_dim1] += tr;
00493     wr[en] = hr[en + en * hr_dim1];
00494     hi[en + en * hi_dim1] += ti;
00495     wi[en] = hi[en + en * hi_dim1];
00496     en = enm1;
00497     goto L220;
00498 /*     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND */
00499 /*                VECTORS OF UPPER TRIANGULAR FORM .......... */
00500 L680:
00501     norm = 0.;
00502 
00503     i__1 = *n;
00504     for (i__ = 1; i__ <= i__1; ++i__) {
00505 
00506         i__2 = *n;
00507         for (j = i__; j <= i__2; ++j) {
00508             tr = (d__1 = hr[i__ + j * hr_dim1], abs(d__1)) + (d__2 = hi[i__ + 
00509                     j * hi_dim1], abs(d__2));
00510             if (tr > norm) {
00511                 norm = tr;
00512             }
00513 /* L720: */
00514         }
00515     }
00516 
00517     if (*n == 1 || norm == 0.) {
00518         goto L1001;
00519     }
00520 /*     .......... FOR EN=N STEP -1 UNTIL 2 DO -- .......... */
00521     i__2 = *n;
00522     for (nn = 2; nn <= i__2; ++nn) {
00523         en = *n + 2 - nn;
00524         xr = wr[en];
00525         xi = wi[en];
00526         hr[en + en * hr_dim1] = 1.;
00527         hi[en + en * hi_dim1] = 0.;
00528         enm1 = en - 1;
00529 /*     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
00530         i__1 = enm1;
00531         for (ii = 1; ii <= i__1; ++ii) {
00532             i__ = en - ii;
00533             zzr = 0.;
00534             zzi = 0.;
00535             ip1 = i__ + 1;
00536 
00537             i__3 = en;
00538             for (j = ip1; j <= i__3; ++j) {
00539                 zzr = zzr + hr[i__ + j * hr_dim1] * hr[j + en * hr_dim1] - hi[
00540                         i__ + j * hi_dim1] * hi[j + en * hi_dim1];
00541                 zzi = zzi + hr[i__ + j * hr_dim1] * hi[j + en * hi_dim1] + hi[
00542                         i__ + j * hi_dim1] * hr[j + en * hr_dim1];
00543 /* L740: */
00544             }
00545 
00546             yr = xr - wr[i__];
00547             yi = xi - wi[i__];
00548             if (yr != 0. || yi != 0.) {
00549                 goto L765;
00550             }
00551             tst1 = norm;
00552             yr = tst1;
00553 L760:
00554             yr *= .01;
00555             tst2 = norm + yr;
00556             if (tst2 > tst1) {
00557                 goto L760;
00558             }
00559 L765:
00560             cdiv_(&zzr, &zzi, &yr, &yi, &hr[i__ + en * hr_dim1], &hi[i__ + en 
00561                     * hi_dim1]);
00562 /*     .......... OVERFLOW CONTROL .......... */
00563             tr = (d__1 = hr[i__ + en * hr_dim1], abs(d__1)) + (d__2 = hi[i__ 
00564                     + en * hi_dim1], abs(d__2));
00565             if (tr == 0.) {
00566                 goto L780;
00567             }
00568             tst1 = tr;
00569             tst2 = tst1 + 1. / tst1;
00570             if (tst2 > tst1) {
00571                 goto L780;
00572             }
00573             i__3 = en;
00574             for (j = i__; j <= i__3; ++j) {
00575                 hr[j + en * hr_dim1] /= tr;
00576                 hi[j + en * hi_dim1] /= tr;
00577 /* L770: */
00578             }
00579 
00580 L780:
00581             ;
00582         }
00583 
00584 /* L800: */
00585     }
00586 /*     .......... END BACKSUBSTITUTION .......... */
00587     enm1 = *n - 1;
00588 /*     .......... VECTORS OF ISOLATED ROOTS .......... */
00589     i__2 = enm1;
00590     for (i__ = 1; i__ <= i__2; ++i__) {
00591         if (i__ >= *low && i__ <= *igh) {
00592             goto L840;
00593         }
00594         ip1 = i__ + 1;
00595 
00596         i__1 = *n;
00597         for (j = ip1; j <= i__1; ++j) {
00598             zr[i__ + j * zr_dim1] = hr[i__ + j * hr_dim1];
00599             zi[i__ + j * zi_dim1] = hi[i__ + j * hi_dim1];
00600 /* L820: */
00601         }
00602 
00603 L840:
00604         ;
00605     }
00606 /*     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
00607 /*                VECTORS OF ORIGINAL FULL MATRIX. */
00608 /*                FOR J=N STEP -1 UNTIL LOW+1 DO -- .......... */
00609     i__2 = enm1;
00610     for (jj = *low; jj <= i__2; ++jj) {
00611         j = *n + *low - jj;
00612         m = min(j,*igh);
00613 
00614         i__1 = *igh;
00615         for (i__ = *low; i__ <= i__1; ++i__) {
00616             zzr = 0.;
00617             zzi = 0.;
00618 
00619             i__3 = m;
00620             for (k = *low; k <= i__3; ++k) {
00621                 zzr = zzr + zr[i__ + k * zr_dim1] * hr[k + j * hr_dim1] - zi[
00622                         i__ + k * zi_dim1] * hi[k + j * hi_dim1];
00623                 zzi = zzi + zr[i__ + k * zr_dim1] * hi[k + j * hi_dim1] + zi[
00624                         i__ + k * zi_dim1] * hr[k + j * hr_dim1];
00625 /* L860: */
00626             }
00627 
00628             zr[i__ + j * zr_dim1] = zzr;
00629             zi[i__ + j * zi_dim1] = zzi;
00630 /* L880: */
00631         }
00632     }
00633 
00634     goto L1001;
00635 /*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
00636 /*                CONVERGED AFTER 30*N ITERATIONS .......... */
00637 L1000:
00638     *ierr = en;
00639 L1001:
00640     return 0;
00641 } /* comqr2_ */
00642