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