Doxygen Source Code Documentation
eis_qzval.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | qzval_ (integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta, logical *matz, doublereal *z__) |
Function Documentation
|
Definition at line 8 of file eis_qzval.c. References a, a2, abs, d_sign(), and v1. 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 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_ */ |