Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_qzval.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int qzval_(integer *nm, integer *n, doublereal *a,
00009 doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta,
00010 logical *matz, doublereal *z__)
00011 {
00012
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
00017 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00018
00019
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
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
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
00119 epsb = b[*n + b_dim1];
00120 isw = 1;
00121
00122
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
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
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
00186
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
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
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
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
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
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
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
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
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
00360
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 }
00398