00014 {
00015     /* System generated locals */
00016     integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2, 
00017             i__3;
00018     doublereal d__1, d__2, d__3;
00019 
00020     /* Builtin functions */
00021     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00022 
00023     /* Local variables */
00024     static doublereal epsa, epsb;
00025     static integer i__, j, k, l;
00026     static doublereal r__, s, t, anorm, bnorm;
00027     static integer enorn;
00028     static doublereal a1, a2, a3;
00029     static integer k1, k2, l1;
00030     static doublereal u1, u2, u3, v1, v2, v3, a11, a12, a21, a22, a33, a34, 
00031             a43, a44, b11, b12, b22, b33;
00032     static integer na, ld;
00033     static doublereal b34, b44;
00034     static integer en;
00035     static doublereal ep;
00036     static integer ll;
00037     static doublereal sh;
00038     extern doublereal epslon_(doublereal *);
00039     static logical notlas;
00040     static integer km1, lm1;
00041     static doublereal ani, bni;
00042     static integer ish, itn, its, enm2, lor1;
00043 
00044 
00045 
00046 /*     THIS SUBROUTINE IS THE SECOND STEP OF THE QZ ALGORITHM */
00047 /*     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */
00048 /*     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, */
00049 /*     AS MODIFIED IN TECHNICAL NOTE NASA TN D-7305(1973) BY WARD. */
00050 
00051 /*     THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM */
00052 /*     IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM. */
00053 /*     IT REDUCES THE HESSENBERG MATRIX TO QUASI-TRIANGULAR FORM USING */
00054 /*     ORTHOGONAL TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM */
00055 /*     OF THE OTHER MATRIX.  IT IS USUALLY PRECEDED BY  QZHES  AND */
00056 /*     FOLLOWED BY  QZVAL  AND, POSSIBLY,  QZVEC. */
00057 
00058 /*     ON INPUT */
00059 
00060 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00061 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00062 /*          DIMENSION STATEMENT. */
00063 
00064 /*        N IS THE ORDER OF THE MATRICES. */
00065 
00066 /*        A CONTAINS A REAL UPPER HESSENBERG MATRIX. */
00067 
00068 /*        B CONTAINS A REAL UPPER TRIANGULAR MATRIX. */
00069 
00070 /*        EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. */
00071 /*          EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN */
00072 /*          ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF */
00073 /*          ERROR TIMES THE NORM OF ITS MATRIX.  IF THE INPUT EPS1 IS */
00074 /*          POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE */
00075 /*          IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX.  A */
00076 /*          POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, */
00077 /*          BUT LESS ACCURATE RESULTS. */
00078 
00079 /*        MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS 
00080 */
00081 /*          ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */
00082 /*          EIGENVECTORS, AND TO .FALSE. OTHERWISE. */
00083 
00084 /*        Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE */
00085 /*          TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION */
00086 /*          BY  QZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. */
00087 /*          IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. */
00088 
00089 /*     ON OUTPUT */
00090 
00091 /*        A HAS BEEN REDUCED TO QUASI-TRIANGULAR FORM.  THE ELEMENTS */
00092 /*          BELOW THE FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO */
00093 /*          CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO. */
00094 
00095 /*        B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS */
00096 /*          HAVE BEEN ALTERED.  THE LOCATION B(N,1) IS USED TO STORE */
00097 /*          EPS1 TIMES THE NORM OF B FOR LATER USE BY  QZVAL  AND  QZVEC. 
00098 */
00099 
00100 /*        Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS */
00101 /*          (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE.. */
00102 
00103 /*        IERR IS SET TO */
00104 /*          ZERO       FOR NORMAL RETURN, */
00105 /*          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
00106 /*                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
00107 
00108 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00109 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00110 */
00111 
00112 /*     THIS VERSION DATED AUGUST 1983. */
00113 
00114 /*     ------------------------------------------------------------------ 
00115 */
00116 
00117     /* Parameter adjustments */
00118     z_dim1 = *nm;
00119     z_offset = z_dim1 + 1;
00120     z__ -= z_offset;
00121     b_dim1 = *nm;
00122     b_offset = b_dim1 + 1;
00123     b -= b_offset;
00124     a_dim1 = *nm;
00125     a_offset = a_dim1 + 1;
00126     a -= a_offset;
00127 
00128     /* Function Body */
00129     *ierr = 0;
00130 /*     .......... COMPUTE EPSA,EPSB .......... */
00131     anorm = 0.;
00132     bnorm = 0.;
00133 
00134     i__1 = *n;
00135     for (i__ = 1; i__ <= i__1; ++i__) {
00136         ani = 0.;
00137         if (i__ != 1) {
00138             ani = (d__1 = a[i__ + (i__ - 1) * a_dim1], abs(d__1));
00139         }
00140         bni = 0.;
00141 
00142         i__2 = *n;
00143         for (j = i__; j <= i__2; ++j) {
00144             ani += (d__1 = a[i__ + j * a_dim1], abs(d__1));
00145             bni += (d__1 = b[i__ + j * b_dim1], abs(d__1));
00146 /* L20: */
00147         }
00148 
00149         if (ani > anorm) {
00150             anorm = ani;
00151         }
00152         if (bni > bnorm) {
00153             bnorm = bni;
00154         }
00155 /* L30: */
00156     }
00157 
00158     if (anorm == 0.) {
00159         anorm = 1.;
00160     }
00161     if (bnorm == 0.) {
00162         bnorm = 1.;
00163     }
00164     ep = *eps1;
00165     if (ep > 0.) {
00166         goto L50;
00167     }
00168 /*     .......... USE ROUNDOFF LEVEL IF EPS1 IS ZERO .......... */
00169     ep = epslon_(&c_b5);
00170 L50:
00171     epsa = ep * anorm;
00172     epsb = ep * bnorm;
00173 /*     .......... REDUCE A TO QUASI-TRIANGULAR FORM, WHILE */
00174 /*                KEEPING B TRIANGULAR .......... */
00175     lor1 = 1;
00176     enorn = *n;
00177     en = *n;
00178     itn = *n * 30;
00179 /*     .......... BEGIN QZ STEP .......... */
00180 L60:
00181     if (en <= 2) {
00182         goto L1001;
00183     }
00184     if (! (*matz)) {
00185         enorn = en;
00186     }
00187     its = 0;
00188     na = en - 1;
00189     enm2 = na - 1;
00190 L70:
00191     ish = 2;
00192 /*     .......... CHECK FOR CONVERGENCE OR REDUCIBILITY. */
00193 /*                FOR L=EN STEP -1 UNTIL 1 DO -- .......... */
00194     i__1 = en;
00195     for (ll = 1; ll <= i__1; ++ll) {
00196         lm1 = en - ll;
00197         l = lm1 + 1;
00198         if (l == 1) {
00199             goto L95;
00200         }
00201         if ((d__1 = a[l + lm1 * a_dim1], abs(d__1)) <= epsa) {
00202             goto L90;
00203         }
00204 /* L80: */
00205     }
00206 
00207 L90:
00208     a[l + lm1 * a_dim1] = 0.;
00209     if (l < na) {
00210         goto L95;
00211     }
00212 /*     .......... 1-BY-1 OR 2-BY-2 BLOCK ISOLATED .......... */
00213     en = lm1;
00214     goto L60;
00215 /*     .......... CHECK FOR SMALL TOP OF B .......... */
00216 L95:
00217     ld = l;
00218 L100:
00219     l1 = l + 1;
00220     b11 = b[l + l * b_dim1];
00221     if (abs(b11) > epsb) {
00222         goto L120;
00223     }
00224     b[l + l * b_dim1] = 0.;
00225     s = (d__1 = a[l + l * a_dim1], abs(d__1)) + (d__2 = a[l1 + l * a_dim1], 
00226             abs(d__2));
00227     u1 = a[l + l * a_dim1] / s;
00228     u2 = a[l1 + l * a_dim1] / s;
00229     d__1 = sqrt(u1 * u1 + u2 * u2);
00230     r__ = d_sign(&d__1, &u1);
00231     v1 = -(u1 + r__) / r__;
00232     v2 = -u2 / r__;
00233     u2 = v2 / v1;
00234 
00235     i__1 = enorn;
00236     for (j = l; j <= i__1; ++j) {
00237         t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1];
00238         a[l + j * a_dim1] += t * v1;
00239         a[l1 + j * a_dim1] += t * v2;
00240         t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1];
00241         b[l + j * b_dim1] += t * v1;
00242         b[l1 + j * b_dim1] += t * v2;
00243 /* L110: */
00244     }
00245 
00246     if (l != 1) {
00247         a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
00248     }
00249     lm1 = l;
00250     l = l1;
00251     goto L90;
00252 L120:
00253     a11 = a[l + l * a_dim1] / b11;
00254     a21 = a[l1 + l * a_dim1] / b11;
00255     if (ish == 1) {
00256         goto L140;
00257     }
00258 /*     .......... ITERATION STRATEGY .......... */
00259     if (itn == 0) {
00260         goto L1000;
00261     }
00262     if (its == 10) {
00263         goto L155;
00264     }
00265 /*     .......... DETERMINE TYPE OF SHIFT .......... */
00266     b22 = b[l1 + l1 * b_dim1];
00267     if (abs(b22) < epsb) {
00268         b22 = epsb;
00269     }
00270     b33 = b[na + na * b_dim1];
00271     if (abs(b33) < epsb) {
00272         b33 = epsb;
00273     }
00274     b44 = b[en + en * b_dim1];
00275     if (abs(b44) < epsb) {
00276         b44 = epsb;
00277     }
00278     a33 = a[na + na * a_dim1] / b33;
00279     a34 = a[na + en * a_dim1] / b44;
00280     a43 = a[en + na * a_dim1] / b33;
00281     a44 = a[en + en * a_dim1] / b44;
00282     b34 = b[na + en * b_dim1] / b44;
00283     t = (a43 * b34 - a33 - a44) * .5;
00284     r__ = t * t + a34 * a43 - a33 * a44;
00285     if (r__ < 0.) {
00286         goto L150;
00287     }
00288 /*     .......... DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A .......... */
00289     ish = 1;
00290     r__ = sqrt(r__);
00291     sh = -t + r__;
00292     s = -t - r__;
00293     if ((d__1 = s - a44, abs(d__1)) < (d__2 = sh - a44, abs(d__2))) {
00294         sh = s;
00295     }
00296 /*     .......... LOOK FOR TWO CONSECUTIVE SMALL */
00297 /*                SUB-DIAGONAL ELEMENTS OF A. */
00298 /*                FOR L=EN-2 STEP -1 UNTIL LD DO -- .......... */
00299     i__1 = enm2;
00300     for (ll = ld; ll <= i__1; ++ll) {
00301         l = enm2 + ld - ll;
00302         if (l == ld) {
00303             goto L140;
00304         }
00305         lm1 = l - 1;
00306         l1 = l + 1;
00307         t = a[l + l * a_dim1];
00308         if ((d__1 = b[l + l * b_dim1], abs(d__1)) > epsb) {
00309             t -= sh * b[l + l * b_dim1];
00310         }
00311         if ((d__1 = a[l + lm1 * a_dim1], abs(d__1)) <= (d__2 = t / a[l1 + l * 
00312                 a_dim1], abs(d__2)) * epsa) {
00313             goto L100;
00314         }
00315 /* L130: */
00316     }
00317 
00318 L140:
00319     a1 = a11 - sh;
00320     a2 = a21;
00321     if (l != ld) {
00322         a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
00323     }
00324     goto L160;
00325 /*     .......... DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A .......... */
00326 L150:
00327     a12 = a[l + l1 * a_dim1] / b22;
00328     a22 = a[l1 + l1 * a_dim1] / b22;
00329     b12 = b[l + l1 * b_dim1] / b22;
00330     a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) / a21 + 
00331             a12 - a11 * b12;
00332     a2 = a22 - a11 - a21 * b12 - (a33 - a11) - (a44 - a11) + a43 * b34;
00333     a3 = a[l1 + 1 + l1 * a_dim1] / b22;
00334     goto L160;
00335 /*     .......... AD HOC SHIFT .......... */
00336 L155:
00337     a1 = 0.;
00338     a2 = 1.;
00339     a3 = 1.1605;
00340 L160:
00341     ++its;
00342     --itn;
00343     if (! (*matz)) {
00344         lor1 = ld;
00345     }
00346 /*     .......... MAIN LOOP .......... */
00347     i__1 = na;
00348     for (k = l; k <= i__1; ++k) {
00349         notlas = k != na && ish == 2;
00350         k1 = k + 1;
00351         k2 = k + 2;
00352 /* Computing MAX */
00353         i__2 = k - 1;
00354         km1 = max(i__2,l);
00355 /* Computing MIN */
00356         i__2 = en, i__3 = k1 + ish;
00357         ll = min(i__2,i__3);
00358         if (notlas) {
00359             goto L190;
00360         }
00361 /*     .......... ZERO A(K+1,K-1) .......... */
00362         if (k == l) {
00363             goto L170;
00364         }
00365         a1 = a[k + km1 * a_dim1];
00366         a2 = a[k1 + km1 * a_dim1];
00367 L170:
00368         s = abs(a1) + abs(a2);
00369         if (s == 0.) {
00370             goto L70;
00371         }
00372         u1 = a1 / s;
00373         u2 = a2 / s;
00374         d__1 = sqrt(u1 * u1 + u2 * u2);
00375         r__ = d_sign(&d__1, &u1);
00376         v1 = -(u1 + r__) / r__;
00377         v2 = -u2 / r__;
00378         u2 = v2 / v1;
00379 
00380         i__2 = enorn;
00381         for (j = km1; j <= i__2; ++j) {
00382             t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1];
00383             a[k + j * a_dim1] += t * v1;
00384             a[k1 + j * a_dim1] += t * v2;
00385             t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1];
00386             b[k + j * b_dim1] += t * v1;
00387             b[k1 + j * b_dim1] += t * v2;
00388 /* L180: */
00389         }
00390 
00391         if (k != l) {
00392             a[k1 + km1 * a_dim1] = 0.;
00393         }
00394         goto L240;
00395 /*     .......... ZERO A(K+1,K-1) AND A(K+2,K-1) .......... */
00396 L190:
00397         if (k == l) {
00398             goto L200;
00399         }
00400         a1 = a[k + km1 * a_dim1];
00401         a2 = a[k1 + km1 * a_dim1];
00402         a3 = a[k2 + km1 * a_dim1];
00403 L200:
00404         s = abs(a1) + abs(a2) + abs(a3);
00405         if (s == 0.) {
00406             goto L260;
00407         }
00408         u1 = a1 / s;
00409         u2 = a2 / s;
00410         u3 = a3 / s;
00411         d__1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
00412         r__ = d_sign(&d__1, &u1);
00413         v1 = -(u1 + r__) / r__;
00414         v2 = -u2 / r__;
00415         v3 = -u3 / r__;
00416         u2 = v2 / v1;
00417         u3 = v3 / v1;
00418 
00419         i__2 = enorn;
00420         for (j = km1; j <= i__2; ++j) {
00421             t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1] + u3 * a[k2 + j * 
00422                     a_dim1];
00423             a[k + j * a_dim1] += t * v1;
00424             a[k1 + j * a_dim1] += t * v2;
00425             a[k2 + j * a_dim1] += t * v3;
00426             t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1] + u3 * b[k2 + j * 
00427                     b_dim1];
00428             b[k + j * b_dim1] += t * v1;
00429             b[k1 + j * b_dim1] += t * v2;
00430             b[k2 + j * b_dim1] += t * v3;
00431 /* L210: */
00432         }
00433 
00434         if (k == l) {
00435             goto L220;
00436         }
00437         a[k1 + km1 * a_dim1] = 0.;
00438         a[k2 + km1 * a_dim1] = 0.;
00439 /*     .......... ZERO B(K+2,K+1) AND B(K+2,K) .......... */
00440 L220:
00441         s = (d__1 = b[k2 + k2 * b_dim1], abs(d__1)) + (d__2 = b[k2 + k1 * 
00442                 b_dim1], abs(d__2)) + (d__3 = b[k2 + k * b_dim1], abs(d__3));
00443         if (s == 0.) {
00444             goto L240;
00445         }
00446         u1 = b[k2 + k2 * b_dim1] / s;
00447         u2 = b[k2 + k1 * b_dim1] / s;
00448         u3 = b[k2 + k * b_dim1] / s;
00449         d__1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
00450         r__ = d_sign(&d__1, &u1);
00451         v1 = -(u1 + r__) / r__;
00452         v2 = -u2 / r__;
00453         v3 = -u3 / r__;
00454         u2 = v2 / v1;
00455         u3 = v3 / v1;
00456 
00457         i__2 = ll;
00458         for (i__ = lor1; i__ <= i__2; ++i__) {
00459             t = a[i__ + k2 * a_dim1] + u2 * a[i__ + k1 * a_dim1] + u3 * a[i__ 
00460                     + k * a_dim1];
00461             a[i__ + k2 * a_dim1] += t * v1;
00462             a[i__ + k1 * a_dim1] += t * v2;
00463             a[i__ + k * a_dim1] += t * v3;
00464             t = b[i__ + k2 * b_dim1] + u2 * b[i__ + k1 * b_dim1] + u3 * b[i__ 
00465                     + k * b_dim1];
00466             b[i__ + k2 * b_dim1] += t * v1;
00467             b[i__ + k1 * b_dim1] += t * v2;
00468             b[i__ + k * b_dim1] += t * v3;
00469 /* L230: */
00470         }
00471 
00472         b[k2 + k * b_dim1] = 0.;
00473         b[k2 + k1 * b_dim1] = 0.;
00474         if (! (*matz)) {
00475             goto L240;
00476         }
00477 
00478         i__2 = *n;
00479         for (i__ = 1; i__ <= i__2; ++i__) {
00480             t = z__[i__ + k2 * z_dim1] + u2 * z__[i__ + k1 * z_dim1] + u3 * 
00481                     z__[i__ + k * z_dim1];
00482             z__[i__ + k2 * z_dim1] += t * v1;
00483             z__[i__ + k1 * z_dim1] += t * v2;
00484             z__[i__ + k * z_dim1] += t * v3;
00485 /* L235: */
00486         }
00487 /*     .......... ZERO B(K+1,K) .......... */
00488 L240:
00489         s = (d__1 = b[k1 + k1 * b_dim1], abs(d__1)) + (d__2 = b[k1 + k * 
00490                 b_dim1], abs(d__2));
00491         if (s == 0.) {
00492             goto L260;
00493         }
00494         u1 = b[k1 + k1 * b_dim1] / s;
00495         u2 = b[k1 + k * b_dim1] / s;
00496         d__1 = sqrt(u1 * u1 + u2 * u2);
00497         r__ = d_sign(&d__1, &u1);
00498         v1 = -(u1 + r__) / r__;
00499         v2 = -u2 / r__;
00500         u2 = v2 / v1;
00501 
00502         i__2 = ll;
00503         for (i__ = lor1; i__ <= i__2; ++i__) {
00504             t = a[i__ + k1 * a_dim1] + u2 * a[i__ + k * a_dim1];
00505             a[i__ + k1 * a_dim1] += t * v1;
00506             a[i__ + k * a_dim1] += t * v2;
00507             t = b[i__ + k1 * b_dim1] + u2 * b[i__ + k * b_dim1];
00508             b[i__ + k1 * b_dim1] += t * v1;
00509             b[i__ + k * b_dim1] += t * v2;
00510 /* L250: */
00511         }
00512 
00513         b[k1 + k * b_dim1] = 0.;
00514         if (! (*matz)) {
00515             goto L260;
00516         }
00517 
00518         i__2 = *n;
00519         for (i__ = 1; i__ <= i__2; ++i__) {
00520             t = z__[i__ + k1 * z_dim1] + u2 * z__[i__ + k * z_dim1];
00521             z__[i__ + k1 * z_dim1] += t * v1;
00522             z__[i__ + k * z_dim1] += t * v2;
00523 /* L255: */
00524         }
00525 
00526 L260:
00527         ;
00528     }
00529 /*     .......... END QZ STEP .......... */
00530     goto L70;
00531 /*     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
00532 /*                CONVERGED AFTER 30*N ITERATIONS .......... */
00533 L1000:
00534     *ierr = en;
00535 /*     .......... SAVE EPSB FOR USE BY QZVAL AND QZVEC .......... */
00536 L1001:
00537     if (*n > 1) {
00538         b[*n + b_dim1] = epsb;
00539     }
00540     return 0;
00541 } /* qzit_ */