Doxygen Source Code Documentation
eis_qzvec.c File Reference
#include "f2c.h"Go to the source code of this file.
Functions | |
| int | qzvec_ (integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta, doublereal *z__) |
Function Documentation
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 8 of file eis_qzvec.c. Referenced by rgg_().
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_ */
|
