Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_qzvec.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int qzvec_(integer *nm, integer *n, doublereal *a,
00009 doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta,
00010 doublereal *z__)
00011 {
00012
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
00018 double sqrt(doublereal);
00019
00020
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
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 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
00113 epsb = b[*n + b_dim1];
00114 isw = 1;
00115
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
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
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
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
00163 t = w;
00164 if (w == 0.) {
00165 t = epsb;
00166 }
00167 b[i__ + en * b_dim1] = -r__ / t;
00168 goto L700;
00169
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
00190 goto L800;
00191
00192 L710:
00193 m = na;
00194 almr = alfr[m];
00195 almi = alfi[m];
00196 betm = beta[m];
00197
00198
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
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
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
00245 tr = -ra;
00246 ti = -sa;
00247 L773:
00248 dr = w;
00249 di = w1;
00250
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
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
00311 L795:
00312 isw = 3 - isw;
00313 L800:
00314 ;
00315 }
00316
00317
00318
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
00330 zz += z__[i__ + k * z_dim1] * b[k + j * b_dim1];
00331 }
00332
00333 z__[i__ + j * z_dim1] = zz;
00334
00335 }
00336 }
00337
00338
00339
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
00356 }
00357
00358 i__1 = *n;
00359 for (i__ = 1; i__ <= i__1; ++i__) {
00360
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
00373 d__1 = z__[i__ + (j - 1) * z_dim1] / r__;
00374
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
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
00389 }
00390
00391 L945:
00392 isw = 3 - isw;
00393 L950:
00394 ;
00395 }
00396
00397 return 0;
00398 }
00399