00001 /* qzvec.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 qzvec_(integer *nm, integer *n, doublereal *a, 
00009         doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta, 
00010         doublereal *z__)
00011 {
00012     /* System generated locals */
00013     integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2, 
00014             i__3;
00015     doublereal d__1, d__2;
00016 
00017     /* Builtin functions */
00018     double sqrt(doublereal);
00019 
00020     /* Local variables */
00021     static doublereal alfm, almi, betm, epsb, almr, d__;
00022     static integer i__, j, k, m;
00023     static doublereal q, r__, s, t, w, x, y, t1, t2, w1, x1, z1, di;
00024     static integer na, ii, en, jj;
00025     static doublereal ra, dr, sa;
00026     static integer nn;
00027     static doublereal ti, rr, tr, zz;
00028     static integer isw, enm2;
00029 
00030 
00031 
00032 /*     THIS SUBROUTINE IS THE OPTIONAL FOURTH STEP OF THE QZ ALGORITHM */
00033 /*     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */
00034 /*     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */
00035 
00036 /*     THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM IN */
00037 /*     QUASI-TRIANGULAR FORM (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS TO */
00038 /*     A PAIR OF COMPLEX EIGENVALUES) AND THE OTHER IN UPPER TRIANGULAR */
00039 /*     FORM.  IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM AND 
00040 */
00041 /*     TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. */
00042 /*     IT IS USUALLY PRECEDED BY  QZHES,  QZIT, AND  QZVAL. */
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 MATRICES. */
00051 
00052 /*        A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. */
00053 
00054 /*        B CONTAINS A REAL UPPER TRIANGULAR MATRIX.  IN ADDITION, */
00055 /*          LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) */
00056 /*          COMPUTED AND SAVED IN  QZIT. */
00057 
00058 /*        ALFR, ALFI, AND BETA  ARE VECTORS WITH COMPONENTS WHOSE */
00059 /*          RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED */
00060 /*          EIGENVALUES.  THEY ARE USUALLY OBTAINED FROM  QZVAL. */
00061 
00062 /*        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */
00063 /*          REDUCTIONS BY  QZHES,  QZIT, AND  QZVAL, IF PERFORMED. */
00064 /*          IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE */
00065 /*          DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX. */
00066 
00067 /*     ON OUTPUT */
00068 
00069 /*        A IS UNALTERED.  ITS SUBDIAGONAL ELEMENTS PROVIDE INFORMATION */
00070 /*           ABOUT THE STORAGE OF THE COMPLEX EIGENVECTORS. */
00071 
00072 /*        B HAS BEEN DESTROYED. */
00073 
00074 /*        ALFR, ALFI, AND BETA ARE UNALTERED. */
00075 
00076 /*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. */
00077 /*          IF ALFI(I) .EQ. 0.0, THE I-TH EIGENVALUE IS REAL AND */
00078 /*            THE I-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. */
00079 /*          IF ALFI(I) .NE. 0.0, THE I-TH EIGENVALUE IS COMPLEX. */
00080 /*            IF ALFI(I) .GT. 0.0, THE EIGENVALUE IS THE FIRST OF */
00081 /*              A COMPLEX PAIR AND THE I-TH AND (I+1)-TH COLUMNS */
00082 /*              OF Z CONTAIN ITS EIGENVECTOR. */
00083 /*            IF ALFI(I) .LT. 0.0, THE EIGENVALUE IS THE SECOND OF */
00084 /*              A COMPLEX PAIR AND THE (I-1)-TH AND I-TH COLUMNS */
00085 /*              OF Z CONTAIN THE CONJUGATE OF ITS EIGENVECTOR. */
00086 /*          EACH EIGENVECTOR IS NORMALIZED SO THAT THE MODULUS */
00087 /*          OF ITS LARGEST COMPONENT IS 1.0 . */
00088 
00089 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00090 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00091 */
00092 
00093 /*     THIS VERSION DATED AUGUST 1983. */
00094 
00095 /*     ------------------------------------------------------------------ 
00096 */
00097 
00098     /* Parameter adjustments */
00099     z_dim1 = *nm;
00100     z_offset = z_dim1 + 1;
00101     z__ -= z_offset;
00102     --beta;
00103     --alfi;
00104     --alfr;
00105     b_dim1 = *nm;
00106     b_offset = b_dim1 + 1;
00107     b -= b_offset;
00108     a_dim1 = *nm;
00109     a_offset = a_dim1 + 1;
00110     a -= a_offset;
00111 
00112     /* Function Body */
00113     epsb = b[*n + b_dim1];
00114     isw = 1;
00115 /*     .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... */
00116     i__1 = *n;
00117     for (nn = 1; nn <= i__1; ++nn) {
00118         en = *n + 1 - nn;
00119         na = en - 1;
00120         if (isw == 2) {
00121             goto L795;
00122         }
00123         if (alfi[en] != 0.) {
00124             goto L710;
00125         }
00126 /*     .......... REAL VECTOR .......... */
00127         m = en;
00128         b[en + en * b_dim1] = 1.;
00129         if (na == 0) {
00130             goto L800;
00131         }
00132         alfm = alfr[m];
00133         betm = beta[m];
00134 /*     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
00135         i__2 = na;
00136         for (ii = 1; ii <= i__2; ++ii) {
00137             i__ = en - ii;
00138             w = betm * a[i__ + i__ * a_dim1] - alfm * b[i__ + i__ * b_dim1];
00139             r__ = 0.;
00140 
00141             i__3 = en;
00142             for (j = m; j <= i__3; ++j) {
00143 /* L610: */
00144                 r__ += (betm * a[i__ + j * a_dim1] - alfm * b[i__ + j * 
00145                         b_dim1]) * b[j + en * b_dim1];
00146             }
00147 
00148             if (i__ == 1 || isw == 2) {
00149                 goto L630;
00150             }
00151             if (betm * a[i__ + (i__ - 1) * a_dim1] == 0.) {
00152                 goto L630;
00153             }
00154             zz = w;
00155             s = r__;
00156             goto L690;
00157 L630:
00158             m = i__;
00159             if (isw == 2) {
00160                 goto L640;
00161             }
00162 /*     .......... REAL 1-BY-1 BLOCK .......... */
00163             t = w;
00164             if (w == 0.) {
00165                 t = epsb;
00166             }
00167             b[i__ + en * b_dim1] = -r__ / t;
00168             goto L700;
00169 /*     .......... REAL 2-BY-2 BLOCK .......... */
00170 L640:
00171             x = betm * a[i__ + (i__ + 1) * a_dim1] - alfm * b[i__ + (i__ + 1) 
00172                     * b_dim1];
00173             y = betm * a[i__ + 1 + i__ * a_dim1];
00174             q = w * zz - x * y;
00175             t = (x * s - zz * r__) / q;
00176             b[i__ + en * b_dim1] = t;
00177             if (abs(x) <= abs(zz)) {
00178                 goto L650;
00179             }
00180             b[i__ + 1 + en * b_dim1] = (-r__ - w * t) / x;
00181             goto L690;
00182 L650:
00183             b[i__ + 1 + en * b_dim1] = (-s - y * t) / zz;
00184 L690:
00185             isw = 3 - isw;
00186 L700:
00187             ;
00188         }
00189 /*     .......... END REAL VECTOR .......... */
00190         goto L800;
00191 /*     .......... COMPLEX VECTOR .......... */
00192 L710:
00193         m = na;
00194         almr = alfr[m];
00195         almi = alfi[m];
00196         betm = beta[m];
00197 /*     .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT */
00198 /*                EIGENVECTOR MATRIX IS TRIANGULAR .......... */
00199         y = betm * a[en + na * a_dim1];
00200         b[na + na * b_dim1] = -almi * b[en + en * b_dim1] / y;
00201         b[na + en * b_dim1] = (almr * b[en + en * b_dim1] - betm * a[en + en *
00202                  a_dim1]) / y;
00203         b[en + na * b_dim1] = 0.;
00204         b[en + en * b_dim1] = 1.;
00205         enm2 = na - 1;
00206         if (enm2 == 0) {
00207             goto L795;
00208         }
00209 /*     .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... */
00210         i__2 = enm2;
00211         for (ii = 1; ii <= i__2; ++ii) {
00212             i__ = na - ii;
00213             w = betm * a[i__ + i__ * a_dim1] - almr * b[i__ + i__ * b_dim1];
00214             w1 = -almi * b[i__ + i__ * b_dim1];
00215             ra = 0.;
00216             sa = 0.;
00217 
00218             i__3 = en;
00219             for (j = m; j <= i__3; ++j) {
00220                 x = betm * a[i__ + j * a_dim1] - almr * b[i__ + j * b_dim1];
00221                 x1 = -almi * b[i__ + j * b_dim1];
00222                 ra = ra + x * b[j + na * b_dim1] - x1 * b[j + en * b_dim1];
00223                 sa = sa + x * b[j + en * b_dim1] + x1 * b[j + na * b_dim1];
00224 /* L760: */
00225             }
00226 
00227             if (i__ == 1 || isw == 2) {
00228                 goto L770;
00229             }
00230             if (betm * a[i__ + (i__ - 1) * a_dim1] == 0.) {
00231                 goto L770;
00232             }
00233             zz = w;
00234             z1 = w1;
00235             r__ = ra;
00236             s = sa;
00237             isw = 2;
00238             goto L790;
00239 L770:
00240             m = i__;
00241             if (isw == 2) {
00242                 goto L780;
00243             }
00244 /*     .......... COMPLEX 1-BY-1 BLOCK .......... */
00245             tr = -ra;
00246             ti = -sa;
00247 L773:
00248             dr = w;
00249             di = w1;
00250 /*     .......... COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR,DI) .....
00251 ..... */
00252 L775:
00253             if (abs(di) > abs(dr)) {
00254                 goto L777;
00255             }
00256             rr = di / dr;
00257             d__ = dr + di * rr;
00258             t1 = (tr + ti * rr) / d__;
00259             t2 = (ti - tr * rr) / d__;
00260             switch (isw) {
00261                 case 1:  goto L787;
00262                 case 2:  goto L782;
00263             }
00264 L777:
00265             rr = dr / di;
00266             d__ = dr * rr + di;
00267             t1 = (tr * rr + ti) / d__;
00268             t2 = (ti * rr - tr) / d__;
00269             switch (isw) {
00270                 case 1:  goto L787;
00271                 case 2:  goto L782;
00272             }
00273 /*     .......... COMPLEX 2-BY-2 BLOCK .......... */
00274 L780:
00275             x = betm * a[i__ + (i__ + 1) * a_dim1] - almr * b[i__ + (i__ + 1) 
00276                     * b_dim1];
00277             x1 = -almi * b[i__ + (i__ + 1) * b_dim1];
00278             y = betm * a[i__ + 1 + i__ * a_dim1];
00279             tr = y * ra - w * r__ + w1 * s;
00280             ti = y * sa - w * s - w1 * r__;
00281             dr = w * zz - w1 * z1 - x * y;
00282             di = w * z1 + w1 * zz - x1 * y;
00283             if (dr == 0. && di == 0.) {
00284                 dr = epsb;
00285             }
00286             goto L775;
00287 L782:
00288             b[i__ + 1 + na * b_dim1] = t1;
00289             b[i__ + 1 + en * b_dim1] = t2;
00290             isw = 1;
00291             if (abs(y) > abs(w) + abs(w1)) {
00292                 goto L785;
00293             }
00294             tr = -ra - x * b[i__ + 1 + na * b_dim1] + x1 * b[i__ + 1 + en * 
00295                     b_dim1];
00296             ti = -sa - x * b[i__ + 1 + en * b_dim1] - x1 * b[i__ + 1 + na * 
00297                     b_dim1];
00298             goto L773;
00299 L785:
00300             t1 = (-r__ - zz * b[i__ + 1 + na * b_dim1] + z1 * b[i__ + 1 + en *
00301                      b_dim1]) / y;
00302             t2 = (-s - zz * b[i__ + 1 + en * b_dim1] - z1 * b[i__ + 1 + na * 
00303                     b_dim1]) / y;
00304 L787:
00305             b[i__ + na * b_dim1] = t1;
00306             b[i__ + en * b_dim1] = t2;
00307 L790:
00308             ;
00309         }
00310 /*     .......... END COMPLEX VECTOR .......... */
00311 L795:
00312         isw = 3 - isw;
00313 L800:
00314         ;
00315     }
00316 /*     .......... END BACK SUBSTITUTION. */
00317 /*                TRANSFORM TO ORIGINAL COORDINATE SYSTEM. */
00318 /*                FOR J=N STEP -1 UNTIL 1 DO -- .......... */
00319     i__1 = *n;
00320     for (jj = 1; jj <= i__1; ++jj) {
00321         j = *n + 1 - jj;
00322 
00323         i__2 = *n;
00324         for (i__ = 1; i__ <= i__2; ++i__) {
00325             zz = 0.;
00326 
00327             i__3 = j;
00328             for (k = 1; k <= i__3; ++k) {
00329 /* L860: */
00330                 zz += z__[i__ + k * z_dim1] * b[k + j * b_dim1];
00331             }
00332 
00333             z__[i__ + j * z_dim1] = zz;
00334 /* L880: */
00335         }
00336     }
00337 /*     .......... NORMALIZE SO THAT MODULUS OF LARGEST */
00338 /*                COMPONENT OF EACH VECTOR IS 1. */
00339 /*                (ISW IS 1 INITIALLY FROM BEFORE) .......... */
00340     i__2 = *n;
00341     for (j = 1; j <= i__2; ++j) {
00342         d__ = 0.;
00343         if (isw == 2) {
00344             goto L920;
00345         }
00346         if (alfi[j] != 0.) {
00347             goto L945;
00348         }
00349 
00350         i__1 = *n;
00351         for (i__ = 1; i__ <= i__1; ++i__) {
00352             if ((d__1 = z__[i__ + j * z_dim1], abs(d__1)) > d__) {
00353                 d__ = (d__2 = z__[i__ + j * z_dim1], abs(d__2));
00354             }
00355 /* L890: */
00356         }
00357 
00358         i__1 = *n;
00359         for (i__ = 1; i__ <= i__1; ++i__) {
00360 /* L900: */
00361             z__[i__ + j * z_dim1] /= d__;
00362         }
00363 
00364         goto L950;
00365 
00366 L920:
00367         i__1 = *n;
00368         for (i__ = 1; i__ <= i__1; ++i__) {
00369             r__ = (d__1 = z__[i__ + (j - 1) * z_dim1], abs(d__1)) + (d__2 = 
00370                     z__[i__ + j * z_dim1], abs(d__2));
00371             if (r__ != 0.) {
00372 /* Computing 2nd power */
00373                 d__1 = z__[i__ + (j - 1) * z_dim1] / r__;
00374 /* Computing 2nd power */
00375                 d__2 = z__[i__ + j * z_dim1] / r__;
00376                 r__ *= sqrt(d__1 * d__1 + d__2 * d__2);
00377             }
00378             if (r__ > d__) {
00379                 d__ = r__;
00380             }
00381 /* L930: */
00382         }
00383 
00384         i__1 = *n;
00385         for (i__ = 1; i__ <= i__1; ++i__) {
00386             z__[i__ + (j - 1) * z_dim1] /= d__;
00387             z__[i__ + j * z_dim1] /= d__;
00388 /* L940: */
00389         }
00390 
00391 L945:
00392         isw = 3 - isw;
00393 L950:
00394         ;
00395     }
00396 
00397     return 0;
00398 } /* qzvec_ */
00399