00011 {
00012     /* System generated locals */
00013     integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2;
00014     doublereal d__1, d__2, d__3, d__4;
00015 
00016     /* Builtin functions */
00017     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00018 
00019     /* Local variables */
00020     static doublereal epsb, c__, d__, e;
00021     static integer i__, j;
00022     static doublereal r__, s, t, a1, a2, u1, u2, v1, v2, a11, a12, a21, a22, 
00023             b11, b12, b22, di, ei;
00024     static integer na;
00025     static doublereal an, bn;
00026     static integer en;
00027     static doublereal cq, dr;
00028     static integer nn;
00029     static doublereal cz, ti, tr, a1i, a2i, a11i, a12i, a22i, a11r, a12r, 
00030             a22r, sqi, ssi;
00031     static integer isw;
00032     static doublereal sqr, szi, ssr, szr;
00033 
00034 
00035 
00036 /*     THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM */
00037 /*     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */
00038 /*     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */
00039 
00040 /*     THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM */
00041 /*     IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM. */
00042 /*     IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY */
00043 /*     REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX */
00044 /*     EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE */
00045 /*     GENERALIZED EIGENVALUES.  IT IS USUALLY PRECEDED BY  QZHES */
00046 /*     AND  QZIT  AND MAY BE FOLLOWED BY  QZVEC. */
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 MATRICES. */
00055 
00056 /*        A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. */
00057 
00058 /*        B CONTAINS A REAL UPPER TRIANGULAR MATRIX.  IN ADDITION, */
00059 /*          LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) */
00060 /*          COMPUTED AND SAVED IN  QZIT. */
00061 
00062 /*        MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS 
00063 */
00064 /*          ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */
00065 /*          EIGENVECTORS, AND TO .FALSE. OTHERWISE. */
00066 
00067 /*        Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE */
00068 /*          TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES */
00069 /*          AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. */
00070 /*          IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. */
00071 
00072 /*     ON OUTPUT */
00073 
00074 /*        A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX */
00075 /*          IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO */
00076 /*          PAIRS OF COMPLEX EIGENVALUES. */
00077 
00078 /*        B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS */
00079 /*          HAVE BEEN ALTERED.  B(N,1) IS UNALTERED. */
00080 
00081 /*        ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE */
00082 /*          DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE */
00083 /*          OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM */
00084 /*          BY UNITARY TRANSFORMATIONS.  NON-ZERO VALUES OF ALFI OCCUR */
00085 /*          IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE. 
00086 */
00087 
00088 /*        BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B, */
00089 /*          NORMALIZED TO BE REAL AND NON-NEGATIVE.  THE GENERALIZED */
00090 /*          EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA). */
00091 
00092 /*        Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS */
00093 /*          (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE. */
00094 
00095 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00096 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00097 */
00098 
00099 /*     THIS VERSION DATED AUGUST 1983. */
00100 
00101 /*     ------------------------------------------------------------------ 
00102 */
00103 
00104     /* Parameter adjustments */
00105     z_dim1 = *nm;
00106     z_offset = z_dim1 + 1;
00107     z__ -= z_offset;
00108     --beta;
00109     --alfi;
00110     --alfr;
00111     b_dim1 = *nm;
00112     b_offset = b_dim1 + 1;
00113     b -= b_offset;
00114     a_dim1 = *nm;
00115     a_offset = a_dim1 + 1;
00116     a -= a_offset;
00117 
00118     /* Function Body */
00119     epsb = b[*n + b_dim1];
00120     isw = 1;
00121 /*     .......... FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES. */
00122 /*                FOR EN=N STEP -1 UNTIL 1 DO -- .......... */
00123     i__1 = *n;
00124     for (nn = 1; nn <= i__1; ++nn) {
00125         en = *n + 1 - nn;
00126         na = en - 1;
00127         if (isw == 2) {
00128             goto L505;
00129         }
00130         if (en == 1) {
00131             goto L410;
00132         }
00133         if (a[en + na * a_dim1] != 0.) {
00134             goto L420;
00135         }
00136 /*     .......... 1-BY-1 BLOCK, ONE REAL ROOT .......... */
00137 L410:
00138         alfr[en] = a[en + en * a_dim1];
00139         if (b[en + en * b_dim1] < 0.) {
00140             alfr[en] = -alfr[en];
00141         }
00142         beta[en] = (d__1 = b[en + en * b_dim1], abs(d__1));
00143         alfi[en] = 0.;
00144         goto L510;
00145 /*     .......... 2-BY-2 BLOCK .......... */
00146 L420:
00147         if ((d__1 = b[na + na * b_dim1], abs(d__1)) <= epsb) {
00148             goto L455;
00149         }
00150         if ((d__1 = b[en + en * b_dim1], abs(d__1)) > epsb) {
00151             goto L430;
00152         }
00153         a1 = a[en + en * a_dim1];
00154         a2 = a[en + na * a_dim1];
00155         bn = 0.;
00156         goto L435;
00157 L430:
00158         an = (d__1 = a[na + na * a_dim1], abs(d__1)) + (d__2 = a[na + en * 
00159                 a_dim1], abs(d__2)) + (d__3 = a[en + na * a_dim1], abs(d__3)) 
00160                 + (d__4 = a[en + en * a_dim1], abs(d__4));
00161         bn = (d__1 = b[na + na * b_dim1], abs(d__1)) + (d__2 = b[na + en * 
00162                 b_dim1], abs(d__2)) + (d__3 = b[en + en * b_dim1], abs(d__3));
00163         a11 = a[na + na * a_dim1] / an;
00164         a12 = a[na + en * a_dim1] / an;
00165         a21 = a[en + na * a_dim1] / an;
00166         a22 = a[en + en * a_dim1] / an;
00167         b11 = b[na + na * b_dim1] / bn;
00168         b12 = b[na + en * b_dim1] / bn;
00169         b22 = b[en + en * b_dim1] / bn;
00170         e = a11 / b11;
00171         ei = a22 / b22;
00172         s = a21 / (b11 * b22);
00173         t = (a22 - e * b22) / b22;
00174         if (abs(e) <= abs(ei)) {
00175             goto L431;
00176         }
00177         e = ei;
00178         t = (a11 - e * b11) / b11;
00179 L431:
00180         c__ = (t - s * b12) * .5;
00181         d__ = c__ * c__ + s * (a12 - e * b12);
00182         if (d__ < 0.) {
00183             goto L480;
00184         }
00185 /*     .......... TWO REAL ROOTS. */
00186 /*                ZERO BOTH A(EN,NA) AND B(EN,NA) .......... */
00187         d__1 = sqrt(d__);
00188         e += c__ + d_sign(&d__1, &c__);
00189         a11 -= e * b11;
00190         a12 -= e * b12;
00191         a22 -= e * b22;
00192         if (abs(a11) + abs(a12) < abs(a21) + abs(a22)) {
00193             goto L432;
00194         }
00195         a1 = a12;
00196         a2 = a11;
00197         goto L435;
00198 L432:
00199         a1 = a22;
00200         a2 = a21;
00201 /*     .......... CHOOSE AND APPLY REAL Z .......... */
00202 L435:
00203         s = abs(a1) + abs(a2);
00204         u1 = a1 / s;
00205         u2 = a2 / s;
00206         d__1 = sqrt(u1 * u1 + u2 * u2);
00207         r__ = d_sign(&d__1, &u1);
00208         v1 = -(u1 + r__) / r__;
00209         v2 = -u2 / r__;
00210         u2 = v2 / v1;
00211 
00212         i__2 = en;
00213         for (i__ = 1; i__ <= i__2; ++i__) {
00214             t = a[i__ + en * a_dim1] + u2 * a[i__ + na * a_dim1];
00215             a[i__ + en * a_dim1] += t * v1;
00216             a[i__ + na * a_dim1] += t * v2;
00217             t = b[i__ + en * b_dim1] + u2 * b[i__ + na * b_dim1];
00218             b[i__ + en * b_dim1] += t * v1;
00219             b[i__ + na * b_dim1] += t * v2;
00220 /* L440: */
00221         }
00222 
00223         if (! (*matz)) {
00224             goto L450;
00225         }
00226 
00227         i__2 = *n;
00228         for (i__ = 1; i__ <= i__2; ++i__) {
00229             t = z__[i__ + en * z_dim1] + u2 * z__[i__ + na * z_dim1];
00230             z__[i__ + en * z_dim1] += t * v1;
00231             z__[i__ + na * z_dim1] += t * v2;
00232 /* L445: */
00233         }
00234 
00235 L450:
00236         if (bn == 0.) {
00237             goto L475;
00238         }
00239         if (an < abs(e) * bn) {
00240             goto L455;
00241         }
00242         a1 = b[na + na * b_dim1];
00243         a2 = b[en + na * b_dim1];
00244         goto L460;
00245 L455:
00246         a1 = a[na + na * a_dim1];
00247         a2 = a[en + na * a_dim1];
00248 /*     .......... CHOOSE AND APPLY REAL Q .......... */
00249 L460:
00250         s = abs(a1) + abs(a2);
00251         if (s == 0.) {
00252             goto L475;
00253         }
00254         u1 = a1 / s;
00255         u2 = a2 / s;
00256         d__1 = sqrt(u1 * u1 + u2 * u2);
00257         r__ = d_sign(&d__1, &u1);
00258         v1 = -(u1 + r__) / r__;
00259         v2 = -u2 / r__;
00260         u2 = v2 / v1;
00261 
00262         i__2 = *n;
00263         for (j = na; j <= i__2; ++j) {
00264             t = a[na + j * a_dim1] + u2 * a[en + j * a_dim1];
00265             a[na + j * a_dim1] += t * v1;
00266             a[en + j * a_dim1] += t * v2;
00267             t = b[na + j * b_dim1] + u2 * b[en + j * b_dim1];
00268             b[na + j * b_dim1] += t * v1;
00269             b[en + j * b_dim1] += t * v2;
00270 /* L470: */
00271         }
00272 
00273 L475:
00274         a[en + na * a_dim1] = 0.;
00275         b[en + na * b_dim1] = 0.;
00276         alfr[na] = a[na + na * a_dim1];
00277         alfr[en] = a[en + en * a_dim1];
00278         if (b[na + na * b_dim1] < 0.) {
00279             alfr[na] = -alfr[na];
00280         }
00281         if (b[en + en * b_dim1] < 0.) {
00282             alfr[en] = -alfr[en];
00283         }
00284         beta[na] = (d__1 = b[na + na * b_dim1], abs(d__1));
00285         beta[en] = (d__1 = b[en + en * b_dim1], abs(d__1));
00286         alfi[en] = 0.;
00287         alfi[na] = 0.;
00288         goto L505;
00289 /*     .......... TWO COMPLEX ROOTS .......... */
00290 L480:
00291         e += c__;
00292         ei = sqrt(-d__);
00293         a11r = a11 - e * b11;
00294         a11i = ei * b11;
00295         a12r = a12 - e * b12;
00296         a12i = ei * b12;
00297         a22r = a22 - e * b22;
00298         a22i = ei * b22;
00299         if (abs(a11r) + abs(a11i) + abs(a12r) + abs(a12i) < abs(a21) + abs(
00300                 a22r) + abs(a22i)) {
00301             goto L482;
00302         }
00303         a1 = a12r;
00304         a1i = a12i;
00305         a2 = -a11r;
00306         a2i = -a11i;
00307         goto L485;
00308 L482:
00309         a1 = a22r;
00310         a1i = a22i;
00311         a2 = -a21;
00312         a2i = 0.;
00313 /*     .......... CHOOSE COMPLEX Z .......... */
00314 L485:
00315         cz = sqrt(a1 * a1 + a1i * a1i);
00316         if (cz == 0.) {
00317             goto L487;
00318         }
00319         szr = (a1 * a2 + a1i * a2i) / cz;
00320         szi = (a1 * a2i - a1i * a2) / cz;
00321         r__ = sqrt(cz * cz + szr * szr + szi * szi);
00322         cz /= r__;
00323         szr /= r__;
00324         szi /= r__;
00325         goto L490;
00326 L487:
00327         szr = 1.;
00328         szi = 0.;
00329 L490:
00330         if (an < (abs(e) + ei) * bn) {
00331             goto L492;
00332         }
00333         a1 = cz * b11 + szr * b12;
00334         a1i = szi * b12;
00335         a2 = szr * b22;
00336         a2i = szi * b22;
00337         goto L495;
00338 L492:
00339         a1 = cz * a11 + szr * a12;
00340         a1i = szi * a12;
00341         a2 = cz * a21 + szr * a22;
00342         a2i = szi * a22;
00343 /*     .......... CHOOSE COMPLEX Q .......... */
00344 L495:
00345         cq = sqrt(a1 * a1 + a1i * a1i);
00346         if (cq == 0.) {
00347             goto L497;
00348         }
00349         sqr = (a1 * a2 + a1i * a2i) / cq;
00350         sqi = (a1 * a2i - a1i * a2) / cq;
00351         r__ = sqrt(cq * cq + sqr * sqr + sqi * sqi);
00352         cq /= r__;
00353         sqr /= r__;
00354         sqi /= r__;
00355         goto L500;
00356 L497:
00357         sqr = 1.;
00358         sqi = 0.;
00359 /*     .......... COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT */
00360 /*                IF TRANSFORMATIONS WERE APPLIED .......... */
00361 L500:
00362         ssr = sqr * szr + sqi * szi;
00363         ssi = sqr * szi - sqi * szr;
00364         i__ = 1;
00365         tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 + ssr * a22;
00366         ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22;
00367         dr = cq * cz * b11 + cq * szr * b12 + ssr * b22;
00368         di = cq * szi * b12 + ssi * b22;
00369         goto L503;
00370 L502:
00371         i__ = 2;
00372         tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 + cq * cz * a22;
00373         ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21;
00374         dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22;
00375         di = -ssi * b11 - sqi * cz * b12;
00376 L503:
00377         t = ti * dr - tr * di;
00378         j = na;
00379         if (t < 0.) {
00380             j = en;
00381         }
00382         r__ = sqrt(dr * dr + di * di);
00383         beta[j] = bn * r__;
00384         alfr[j] = an * (tr * dr + ti * di) / r__;
00385         alfi[j] = an * t / r__;
00386         if (i__ == 1) {
00387             goto L502;
00388         }
00389 L505:
00390         isw = 3 - isw;
00391 L510:
00392         ;
00393     }
00394     b[*n + b_dim1] = epsb;
00395 
00396     return 0;
00397 } /* qzval_ */