00001 /* qzhes.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 qzhes_(integer *nm, integer *n, doublereal *a, 
00009         doublereal *b, logical *matz, doublereal *z__)
00010 {
00011     /* System generated locals */
00012     integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2, 
00013             i__3;
00014     doublereal d__1, d__2;
00015 
00016     /* Builtin functions */
00017     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00018 
00019     /* Local variables */
00020     static integer i__, j, k, l;
00021     static doublereal r__, s, t;
00022     static integer l1;
00023     static doublereal u1, u2, v1, v2;
00024     static integer lb, nk1, nm1, nm2;
00025     static doublereal rho;
00026 
00027 
00028 
00029 /*     THIS SUBROUTINE IS THE FIRST STEP OF THE QZ ALGORITHM */
00030 /*     FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */
00031 /*     SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */
00032 
00033 /*     THIS SUBROUTINE ACCEPTS A PAIR OF REAL GENERAL MATRICES AND */
00034 /*     REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER */
00035 /*     TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS. */
00036 /*     IT IS USUALLY FOLLOWED BY  QZIT,  QZVAL  AND, POSSIBLY,  QZVEC. */
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 MATRICES. */
00045 
00046 /*        A CONTAINS A REAL GENERAL MATRIX. */
00047 
00048 /*        B CONTAINS A REAL GENERAL MATRIX. */
00049 
00050 /*        MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS 
00051 */
00052 /*          ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */
00053 /*          EIGENVECTORS, AND TO .FALSE. OTHERWISE. */
00054 
00055 /*     ON OUTPUT */
00056 
00057 /*        A HAS BEEN REDUCED TO UPPER HESSENBERG FORM.  THE ELEMENTS */
00058 /*          BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO. */
00059 
00060 /*        B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM.  THE ELEMENTS */
00061 /*          BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO. */
00062 
00063 /*        Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF */
00064 /*          MATZ HAS BEEN SET TO .TRUE.  OTHERWISE, Z IS NOT REFERENCED. 
00065 */
00066 
00067 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00068 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00069 */
00070 
00071 /*     THIS VERSION DATED AUGUST 1983. */
00072 
00073 /*     ------------------------------------------------------------------ 
00074 */
00075 
00076 /*     .......... INITIALIZE Z .......... */
00077     /* Parameter adjustments */
00078     z_dim1 = *nm;
00079     z_offset = z_dim1 + 1;
00080     z__ -= z_offset;
00081     b_dim1 = *nm;
00082     b_offset = b_dim1 + 1;
00083     b -= b_offset;
00084     a_dim1 = *nm;
00085     a_offset = a_dim1 + 1;
00086     a -= a_offset;
00087 
00088     /* Function Body */
00089     if (! (*matz)) {
00090         goto L10;
00091     }
00092 
00093     i__1 = *n;
00094     for (j = 1; j <= i__1; ++j) {
00095 
00096         i__2 = *n;
00097         for (i__ = 1; i__ <= i__2; ++i__) {
00098             z__[i__ + j * z_dim1] = 0.;
00099 /* L2: */
00100         }
00101 
00102         z__[j + j * z_dim1] = 1.;
00103 /* L3: */
00104     }
00105 /*     .......... REDUCE B TO UPPER TRIANGULAR FORM .......... */
00106 L10:
00107     if (*n <= 1) {
00108         goto L170;
00109     }
00110     nm1 = *n - 1;
00111 
00112     i__1 = nm1;
00113     for (l = 1; l <= i__1; ++l) {
00114         l1 = l + 1;
00115         s = 0.;
00116 
00117         i__2 = *n;
00118         for (i__ = l1; i__ <= i__2; ++i__) {
00119             s += (d__1 = b[i__ + l * b_dim1], abs(d__1));
00120 /* L20: */
00121         }
00122 
00123         if (s == 0.) {
00124             goto L100;
00125         }
00126         s += (d__1 = b[l + l * b_dim1], abs(d__1));
00127         r__ = 0.;
00128 
00129         i__2 = *n;
00130         for (i__ = l; i__ <= i__2; ++i__) {
00131             b[i__ + l * b_dim1] /= s;
00132 /* Computing 2nd power */
00133             d__1 = b[i__ + l * b_dim1];
00134             r__ += d__1 * d__1;
00135 /* L25: */
00136         }
00137 
00138         d__1 = sqrt(r__);
00139         r__ = d_sign(&d__1, &b[l + l * b_dim1]);
00140         b[l + l * b_dim1] += r__;
00141         rho = r__ * b[l + l * b_dim1];
00142 
00143         i__2 = *n;
00144         for (j = l1; j <= i__2; ++j) {
00145             t = 0.;
00146 
00147             i__3 = *n;
00148             for (i__ = l; i__ <= i__3; ++i__) {
00149                 t += b[i__ + l * b_dim1] * b[i__ + j * b_dim1];
00150 /* L30: */
00151             }
00152 
00153             t = -t / rho;
00154 
00155             i__3 = *n;
00156             for (i__ = l; i__ <= i__3; ++i__) {
00157                 b[i__ + j * b_dim1] += t * b[i__ + l * b_dim1];
00158 /* L40: */
00159             }
00160 
00161 /* L50: */
00162         }
00163 
00164         i__2 = *n;
00165         for (j = 1; j <= i__2; ++j) {
00166             t = 0.;
00167 
00168             i__3 = *n;
00169             for (i__ = l; i__ <= i__3; ++i__) {
00170                 t += b[i__ + l * b_dim1] * a[i__ + j * a_dim1];
00171 /* L60: */
00172             }
00173 
00174             t = -t / rho;
00175 
00176             i__3 = *n;
00177             for (i__ = l; i__ <= i__3; ++i__) {
00178                 a[i__ + j * a_dim1] += t * b[i__ + l * b_dim1];
00179 /* L70: */
00180             }
00181 
00182 /* L80: */
00183         }
00184 
00185         b[l + l * b_dim1] = -s * r__;
00186 
00187         i__2 = *n;
00188         for (i__ = l1; i__ <= i__2; ++i__) {
00189             b[i__ + l * b_dim1] = 0.;
00190 /* L90: */
00191         }
00192 
00193 L100:
00194         ;
00195     }
00196 /*     .......... REDUCE A TO UPPER HESSENBERG FORM, WHILE */
00197 /*                KEEPING B TRIANGULAR .......... */
00198     if (*n == 2) {
00199         goto L170;
00200     }
00201     nm2 = *n - 2;
00202 
00203     i__1 = nm2;
00204     for (k = 1; k <= i__1; ++k) {
00205         nk1 = nm1 - k;
00206 /*     .......... FOR L=N-1 STEP -1 UNTIL K+1 DO -- .......... */
00207         i__2 = nk1;
00208         for (lb = 1; lb <= i__2; ++lb) {
00209             l = *n - lb;
00210             l1 = l + 1;
00211 /*     .......... ZERO A(L+1,K) .......... */
00212             s = (d__1 = a[l + k * a_dim1], abs(d__1)) + (d__2 = a[l1 + k * 
00213                     a_dim1], abs(d__2));
00214             if (s == 0.) {
00215                 goto L150;
00216             }
00217             u1 = a[l + k * a_dim1] / s;
00218             u2 = a[l1 + k * a_dim1] / s;
00219             d__1 = sqrt(u1 * u1 + u2 * u2);
00220             r__ = d_sign(&d__1, &u1);
00221             v1 = -(u1 + r__) / r__;
00222             v2 = -u2 / r__;
00223             u2 = v2 / v1;
00224 
00225             i__3 = *n;
00226             for (j = k; j <= i__3; ++j) {
00227                 t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1];
00228                 a[l + j * a_dim1] += t * v1;
00229                 a[l1 + j * a_dim1] += t * v2;
00230 /* L110: */
00231             }
00232 
00233             a[l1 + k * a_dim1] = 0.;
00234 
00235             i__3 = *n;
00236             for (j = l; j <= i__3; ++j) {
00237                 t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1];
00238                 b[l + j * b_dim1] += t * v1;
00239                 b[l1 + j * b_dim1] += t * v2;
00240 /* L120: */
00241             }
00242 /*     .......... ZERO B(L+1,L) .......... */
00243             s = (d__1 = b[l1 + l1 * b_dim1], abs(d__1)) + (d__2 = b[l1 + l * 
00244                     b_dim1], abs(d__2));
00245             if (s == 0.) {
00246                 goto L150;
00247             }
00248             u1 = b[l1 + l1 * b_dim1] / s;
00249             u2 = b[l1 + l * b_dim1] / s;
00250             d__1 = sqrt(u1 * u1 + u2 * u2);
00251             r__ = d_sign(&d__1, &u1);
00252             v1 = -(u1 + r__) / r__;
00253             v2 = -u2 / r__;
00254             u2 = v2 / v1;
00255 
00256             i__3 = l1;
00257             for (i__ = 1; i__ <= i__3; ++i__) {
00258                 t = b[i__ + l1 * b_dim1] + u2 * b[i__ + l * b_dim1];
00259                 b[i__ + l1 * b_dim1] += t * v1;
00260                 b[i__ + l * b_dim1] += t * v2;
00261 /* L130: */
00262             }
00263 
00264             b[l1 + l * b_dim1] = 0.;
00265 
00266             i__3 = *n;
00267             for (i__ = 1; i__ <= i__3; ++i__) {
00268                 t = a[i__ + l1 * a_dim1] + u2 * a[i__ + l * a_dim1];
00269                 a[i__ + l1 * a_dim1] += t * v1;
00270                 a[i__ + l * a_dim1] += t * v2;
00271 /* L140: */
00272             }
00273 
00274             if (! (*matz)) {
00275                 goto L150;
00276             }
00277 
00278             i__3 = *n;
00279             for (i__ = 1; i__ <= i__3; ++i__) {
00280                 t = z__[i__ + l1 * z_dim1] + u2 * z__[i__ + l * z_dim1];
00281                 z__[i__ + l1 * z_dim1] += t * v1;
00282                 z__[i__ + l * z_dim1] += t * v2;
00283 /* L145: */
00284             }
00285 
00286 L150:
00287             ;
00288         }
00289 
00290 /* L160: */
00291     }
00292 
00293 L170:
00294     return 0;
00295 } /* qzhes_ */
00296