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