00001 /* comlr.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 comlr_(integer *nm, integer *n, integer *low, integer *
00009         igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, 
00010         integer *ierr)
00011 {
00012     /* System generated locals */
00013     integer hr_dim1, hr_offset, hi_dim1, hi_offset, i__1, i__2;
00014     doublereal d__1, d__2, d__3, d__4;
00015 
00016     /* Local variables */
00017     extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
00018             , doublereal *, doublereal *, doublereal *);
00019     static integer i__, j, l, m, en, ll, mm;
00020     static doublereal si, ti, xi, yi, sr, tr, xr, yr;
00021     static integer im1;
00022     extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 
00023             doublereal *, doublereal *);
00024     static integer mp1, itn, its;
00025     static doublereal zzi, zzr;
00026     static integer enm1;
00027     static doublereal tst1, tst2;
00028 
00029 
00030 
00031 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR, */
00032 /*     NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON. */
00033 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). */
00034 
00035 /*     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX */
00036 /*     UPPER HESSENBERG MATRIX BY THE MODIFIED LR METHOD. */
00037 
00038 /*     ON INPUT */
00039 
00040 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00041 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00042 /*          DIMENSION STATEMENT. */
00043 
00044 /*        N IS THE ORDER OF THE MATRIX. */
00045 
00046 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00047 /*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
00048 /*          SET LOW=1, IGH=N. */
00049 
00050 /*        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
00051 /*          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
00052 /*          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE */
00053 /*          MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY  COMHES, */
00054 /*          IF PERFORMED. */
00055 
00056 /*     ON OUTPUT */
00057 
00058 /*        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
00059 /*          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE */
00060 /*          CALLING  COMLR  IF SUBSEQUENT CALCULATION OF */
00061 /*          EIGENVECTORS IS TO BE PERFORMED. */
00062 
00063 /*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
00064 /*          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR */
00065 /*          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 /*     CALLS CDIV FOR COMPLEX DIVISION. */
00074 /*     CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
00075 
00076 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00077 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00078 */
00079 
00080 /*     THIS VERSION DATED AUGUST 1983. */
00081 
00082 /*     ------------------------------------------------------------------ 
00083 */
00084 
00085     /* Parameter adjustments */
00086     --wi;
00087     --wr;
00088     hi_dim1 = *nm;
00089     hi_offset = hi_dim1 + 1;
00090     hi -= hi_offset;
00091     hr_dim1 = *nm;
00092     hr_offset = hr_dim1 + 1;
00093     hr -= hr_offset;
00094 
00095     /* Function Body */
00096     *ierr = 0;
00097 /*     .......... STORE ROOTS ISOLATED BY CBAL .......... */
00098     i__1 = *n;
00099     for (i__ = 1; i__ <= i__1; ++i__) {
00100         if (i__ >= *low && i__ <= *igh) {
00101             goto L200;
00102         }
00103         wr[i__] = hr[i__ + i__ * hr_dim1];
00104         wi[i__] = hi[i__ + i__ * hi_dim1];
00105 L200:
00106         ;
00107     }
00108 
00109     en = *igh;
00110     tr = 0.;
00111     ti = 0.;
00112     itn = *n * 30;
00113 /*     .......... SEARCH FOR NEXT EIGENVALUE .......... */
00114 L220:
00115     if (en < *low) {
00116         goto L1001;
00117     }
00118     its = 0;
00119     enm1 = en - 1;
00120 /*     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
00121 /*                FOR L=EN STEP -1 UNTIL LOW D0 -- .......... */
00122 L240:
00123     i__1 = en;
00124     for (ll = *low; ll <= i__1; ++ll) {
00125         l = en + *low - ll;
00126         if (l == *low) {
00127             goto L300;
00128         }
00129         tst1 = (d__1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[
00130                 l - 1 + (l - 1) * hi_dim1], abs(d__2)) + (d__3 = hr[l + l * 
00131                 hr_dim1], abs(d__3)) + (d__4 = hi[l + l * hi_dim1], abs(d__4))
00132                 ;
00133         tst2 = tst1 + (d__1 = hr[l + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = 
00134                 hi[l + (l - 1) * hi_dim1], abs(d__2));
00135         if (tst2 == tst1) {
00136             goto L300;
00137         }
00138 /* L260: */
00139     }
00140 /*     .......... FORM SHIFT .......... */
00141 L300:
00142     if (l == en) {
00143         goto L660;
00144     }
00145     if (itn == 0) {
00146         goto L1000;
00147     }
00148     if (its == 10 || its == 20) {
00149         goto L320;
00150     }
00151     sr = hr[en + en * hr_dim1];
00152     si = hi[en + en * hi_dim1];
00153     xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en * 
00154             hi_dim1] * hi[en + enm1 * hi_dim1];
00155     xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en * 
00156             hi_dim1] * hr[en + enm1 * hr_dim1];
00157     if (xr == 0. && xi == 0.) {
00158         goto L340;
00159     }
00160     yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
00161     yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
00162 /* Computing 2nd power */
00163     d__2 = yr;
00164 /* Computing 2nd power */
00165     d__3 = yi;
00166     d__1 = d__2 * d__2 - d__3 * d__3 + xr;
00167     d__4 = yr * 2. * yi + xi;
00168     csroot_(&d__1, &d__4, &zzr, &zzi);
00169     if (yr * zzr + yi * zzi >= 0.) {
00170         goto L310;
00171     }
00172     zzr = -zzr;
00173     zzi = -zzi;
00174 L310:
00175     d__1 = yr + zzr;
00176     d__2 = yi + zzi;
00177     cdiv_(&xr, &xi, &d__1, &d__2, &xr, &xi);
00178     sr -= xr;
00179     si -= xi;
00180     goto L340;
00181 /*     .......... FORM EXCEPTIONAL SHIFT .......... */
00182 L320:
00183     sr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hr[enm1 + (en 
00184             - 2) * hr_dim1], abs(d__2));
00185     si = (d__1 = hi[en + enm1 * hi_dim1], abs(d__1)) + (d__2 = hi[enm1 + (en 
00186             - 2) * hi_dim1], abs(d__2));
00187 
00188 L340:
00189     i__1 = en;
00190     for (i__ = *low; i__ <= i__1; ++i__) {
00191         hr[i__ + i__ * hr_dim1] -= sr;
00192         hi[i__ + i__ * hi_dim1] -= si;
00193 /* L360: */
00194     }
00195 
00196     tr += sr;
00197     ti += si;
00198     ++its;
00199     --itn;
00200 /*     .......... LOOK FOR TWO CONSECUTIVE SMALL */
00201 /*                SUB-DIAGONAL ELEMENTS .......... */
00202     xr = (d__1 = hr[enm1 + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[enm1 + 
00203             enm1 * hi_dim1], abs(d__2));
00204     yr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[en + enm1 * 
00205             hi_dim1], abs(d__2));
00206     zzr = (d__1 = hr[en + en * hr_dim1], abs(d__1)) + (d__2 = hi[en + en * 
00207             hi_dim1], abs(d__2));
00208 /*     .......... FOR M=EN-1 STEP -1 UNTIL L DO -- .......... */
00209     i__1 = enm1;
00210     for (mm = l; mm <= i__1; ++mm) {
00211         m = enm1 + l - mm;
00212         if (m == l) {
00213             goto L420;
00214         }
00215         yi = yr;
00216         yr = (d__1 = hr[m + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m + (
00217                 m - 1) * hi_dim1], abs(d__2));
00218         xi = zzr;
00219         zzr = xr;
00220         xr = (d__1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m 
00221                 - 1 + (m - 1) * hi_dim1], abs(d__2));
00222         tst1 = zzr / yi * (zzr + xr + xi);
00223         tst2 = tst1 + yr;
00224         if (tst2 == tst1) {
00225             goto L420;
00226         }
00227 /* L380: */
00228     }
00229 /*     .......... TRIANGULAR DECOMPOSITION H=L*R .......... */
00230 L420:
00231     mp1 = m + 1;
00232 
00233     i__1 = en;
00234     for (i__ = mp1; i__ <= i__1; ++i__) {
00235         im1 = i__ - 1;
00236         xr = hr[im1 + im1 * hr_dim1];
00237         xi = hi[im1 + im1 * hi_dim1];
00238         yr = hr[i__ + im1 * hr_dim1];
00239         yi = hi[i__ + im1 * hi_dim1];
00240         if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
00241             goto L460;
00242         }
00243 /*     .......... INTERCHANGE ROWS OF HR AND HI .......... */
00244         i__2 = en;
00245         for (j = im1; j <= i__2; ++j) {
00246             zzr = hr[im1 + j * hr_dim1];
00247             hr[im1 + j * hr_dim1] = hr[i__ + j * hr_dim1];
00248             hr[i__ + j * hr_dim1] = zzr;
00249             zzi = hi[im1 + j * hi_dim1];
00250             hi[im1 + j * hi_dim1] = hi[i__ + j * hi_dim1];
00251             hi[i__ + j * hi_dim1] = zzi;
00252 /* L440: */
00253         }
00254 
00255         cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
00256         wr[i__] = 1.;
00257         goto L480;
00258 L460:
00259         cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
00260         wr[i__] = -1.;
00261 L480:
00262         hr[i__ + im1 * hr_dim1] = zzr;
00263         hi[i__ + im1 * hi_dim1] = zzi;
00264 
00265         i__2 = en;
00266         for (j = i__; j <= i__2; ++j) {
00267             hr[i__ + j * hr_dim1] = hr[i__ + j * hr_dim1] - zzr * hr[im1 + j *
00268                      hr_dim1] + zzi * hi[im1 + j * hi_dim1];
00269             hi[i__ + j * hi_dim1] = hi[i__ + j * hi_dim1] - zzr * hi[im1 + j *
00270                      hi_dim1] - zzi * hr[im1 + j * hr_dim1];
00271 /* L500: */
00272         }
00273 
00274 /* L520: */
00275     }
00276 /*     .......... COMPOSITION R*L=H .......... */
00277     i__1 = en;
00278     for (j = mp1; j <= i__1; ++j) {
00279         xr = hr[j + (j - 1) * hr_dim1];
00280         xi = hi[j + (j - 1) * hi_dim1];
00281         hr[j + (j - 1) * hr_dim1] = 0.;
00282         hi[j + (j - 1) * hi_dim1] = 0.;
00283 /*     .......... INTERCHANGE COLUMNS OF HR AND HI, */
00284 /*                IF NECESSARY .......... */
00285         if (wr[j] <= 0.) {
00286             goto L580;
00287         }
00288 
00289         i__2 = j;
00290         for (i__ = l; i__ <= i__2; ++i__) {
00291             zzr = hr[i__ + (j - 1) * hr_dim1];
00292             hr[i__ + (j - 1) * hr_dim1] = hr[i__ + j * hr_dim1];
00293             hr[i__ + j * hr_dim1] = zzr;
00294             zzi = hi[i__ + (j - 1) * hi_dim1];
00295             hi[i__ + (j - 1) * hi_dim1] = hi[i__ + j * hi_dim1];
00296             hi[i__ + j * hi_dim1] = zzi;
00297 /* L540: */
00298         }
00299 
00300 L580:
00301         i__2 = j;
00302         for (i__ = l; i__ <= i__2; ++i__) {
00303             hr[i__ + (j - 1) * hr_dim1] = hr[i__ + (j - 1) * hr_dim1] + xr * 
00304                     hr[i__ + j * hr_dim1] - xi * hi[i__ + j * hi_dim1];
00305             hi[i__ + (j - 1) * hi_dim1] = hi[i__ + (j - 1) * hi_dim1] + xr * 
00306                     hi[i__ + j * hi_dim1] + xi * hr[i__ + j * hr_dim1];
00307 /* L600: */
00308         }
00309 
00310 /* L640: */
00311     }
00312 
00313     goto L240;
00314 /*     .......... A ROOT FOUND .......... */
00315 L660:
00316     wr[en] = hr[en + en * hr_dim1] + tr;
00317     wi[en] = hi[en + en * hi_dim1] + ti;
00318     en = enm1;
00319     goto L220;
00320 /*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
00321 /*                CONVERGED AFTER 30*N ITERATIONS .......... */
00322 L1000:
00323     *ierr = en;
00324 L1001:
00325     return 0;
00326 } /* comlr_ */
00327