Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_qzit.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008
00009
00010 static doublereal c_b5 = 1.;
00011
00012 int qzit_(integer *nm, integer *n, doublereal *a, doublereal
00013 *b, doublereal *eps1, logical *matz, doublereal *z__, integer *ierr)
00014 {
00015
00016 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2,
00017 i__3;
00018 doublereal d__1, d__2, d__3;
00019
00020
00021 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00022
00023
00024 static doublereal epsa, epsb;
00025 static integer i__, j, k, l;
00026 static doublereal r__, s, t, anorm, bnorm;
00027 static integer enorn;
00028 static doublereal a1, a2, a3;
00029 static integer k1, k2, l1;
00030 static doublereal u1, u2, u3, v1, v2, v3, a11, a12, a21, a22, a33, a34,
00031 a43, a44, b11, b12, b22, b33;
00032 static integer na, ld;
00033 static doublereal b34, b44;
00034 static integer en;
00035 static doublereal ep;
00036 static integer ll;
00037 static doublereal sh;
00038 extern doublereal epslon_(doublereal *);
00039 static logical notlas;
00040 static integer km1, lm1;
00041 static doublereal ani, bni;
00042 static integer ish, itn, its, enm2, lor1;
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
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118 z_dim1 = *nm;
00119 z_offset = z_dim1 + 1;
00120 z__ -= z_offset;
00121 b_dim1 = *nm;
00122 b_offset = b_dim1 + 1;
00123 b -= b_offset;
00124 a_dim1 = *nm;
00125 a_offset = a_dim1 + 1;
00126 a -= a_offset;
00127
00128
00129 *ierr = 0;
00130
00131 anorm = 0.;
00132 bnorm = 0.;
00133
00134 i__1 = *n;
00135 for (i__ = 1; i__ <= i__1; ++i__) {
00136 ani = 0.;
00137 if (i__ != 1) {
00138 ani = (d__1 = a[i__ + (i__ - 1) * a_dim1], abs(d__1));
00139 }
00140 bni = 0.;
00141
00142 i__2 = *n;
00143 for (j = i__; j <= i__2; ++j) {
00144 ani += (d__1 = a[i__ + j * a_dim1], abs(d__1));
00145 bni += (d__1 = b[i__ + j * b_dim1], abs(d__1));
00146
00147 }
00148
00149 if (ani > anorm) {
00150 anorm = ani;
00151 }
00152 if (bni > bnorm) {
00153 bnorm = bni;
00154 }
00155
00156 }
00157
00158 if (anorm == 0.) {
00159 anorm = 1.;
00160 }
00161 if (bnorm == 0.) {
00162 bnorm = 1.;
00163 }
00164 ep = *eps1;
00165 if (ep > 0.) {
00166 goto L50;
00167 }
00168
00169 ep = epslon_(&c_b5);
00170 L50:
00171 epsa = ep * anorm;
00172 epsb = ep * bnorm;
00173
00174
00175 lor1 = 1;
00176 enorn = *n;
00177 en = *n;
00178 itn = *n * 30;
00179
00180 L60:
00181 if (en <= 2) {
00182 goto L1001;
00183 }
00184 if (! (*matz)) {
00185 enorn = en;
00186 }
00187 its = 0;
00188 na = en - 1;
00189 enm2 = na - 1;
00190 L70:
00191 ish = 2;
00192
00193
00194 i__1 = en;
00195 for (ll = 1; ll <= i__1; ++ll) {
00196 lm1 = en - ll;
00197 l = lm1 + 1;
00198 if (l == 1) {
00199 goto L95;
00200 }
00201 if ((d__1 = a[l + lm1 * a_dim1], abs(d__1)) <= epsa) {
00202 goto L90;
00203 }
00204
00205 }
00206
00207 L90:
00208 a[l + lm1 * a_dim1] = 0.;
00209 if (l < na) {
00210 goto L95;
00211 }
00212
00213 en = lm1;
00214 goto L60;
00215
00216 L95:
00217 ld = l;
00218 L100:
00219 l1 = l + 1;
00220 b11 = b[l + l * b_dim1];
00221 if (abs(b11) > epsb) {
00222 goto L120;
00223 }
00224 b[l + l * b_dim1] = 0.;
00225 s = (d__1 = a[l + l * a_dim1], abs(d__1)) + (d__2 = a[l1 + l * a_dim1],
00226 abs(d__2));
00227 u1 = a[l + l * a_dim1] / s;
00228 u2 = a[l1 + l * a_dim1] / s;
00229 d__1 = sqrt(u1 * u1 + u2 * u2);
00230 r__ = d_sign(&d__1, &u1);
00231 v1 = -(u1 + r__) / r__;
00232 v2 = -u2 / r__;
00233 u2 = v2 / v1;
00234
00235 i__1 = enorn;
00236 for (j = l; j <= i__1; ++j) {
00237 t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1];
00238 a[l + j * a_dim1] += t * v1;
00239 a[l1 + j * a_dim1] += t * v2;
00240 t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1];
00241 b[l + j * b_dim1] += t * v1;
00242 b[l1 + j * b_dim1] += t * v2;
00243
00244 }
00245
00246 if (l != 1) {
00247 a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
00248 }
00249 lm1 = l;
00250 l = l1;
00251 goto L90;
00252 L120:
00253 a11 = a[l + l * a_dim1] / b11;
00254 a21 = a[l1 + l * a_dim1] / b11;
00255 if (ish == 1) {
00256 goto L140;
00257 }
00258
00259 if (itn == 0) {
00260 goto L1000;
00261 }
00262 if (its == 10) {
00263 goto L155;
00264 }
00265
00266 b22 = b[l1 + l1 * b_dim1];
00267 if (abs(b22) < epsb) {
00268 b22 = epsb;
00269 }
00270 b33 = b[na + na * b_dim1];
00271 if (abs(b33) < epsb) {
00272 b33 = epsb;
00273 }
00274 b44 = b[en + en * b_dim1];
00275 if (abs(b44) < epsb) {
00276 b44 = epsb;
00277 }
00278 a33 = a[na + na * a_dim1] / b33;
00279 a34 = a[na + en * a_dim1] / b44;
00280 a43 = a[en + na * a_dim1] / b33;
00281 a44 = a[en + en * a_dim1] / b44;
00282 b34 = b[na + en * b_dim1] / b44;
00283 t = (a43 * b34 - a33 - a44) * .5;
00284 r__ = t * t + a34 * a43 - a33 * a44;
00285 if (r__ < 0.) {
00286 goto L150;
00287 }
00288
00289 ish = 1;
00290 r__ = sqrt(r__);
00291 sh = -t + r__;
00292 s = -t - r__;
00293 if ((d__1 = s - a44, abs(d__1)) < (d__2 = sh - a44, abs(d__2))) {
00294 sh = s;
00295 }
00296
00297
00298
00299 i__1 = enm2;
00300 for (ll = ld; ll <= i__1; ++ll) {
00301 l = enm2 + ld - ll;
00302 if (l == ld) {
00303 goto L140;
00304 }
00305 lm1 = l - 1;
00306 l1 = l + 1;
00307 t = a[l + l * a_dim1];
00308 if ((d__1 = b[l + l * b_dim1], abs(d__1)) > epsb) {
00309 t -= sh * b[l + l * b_dim1];
00310 }
00311 if ((d__1 = a[l + lm1 * a_dim1], abs(d__1)) <= (d__2 = t / a[l1 + l *
00312 a_dim1], abs(d__2)) * epsa) {
00313 goto L100;
00314 }
00315
00316 }
00317
00318 L140:
00319 a1 = a11 - sh;
00320 a2 = a21;
00321 if (l != ld) {
00322 a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
00323 }
00324 goto L160;
00325
00326 L150:
00327 a12 = a[l + l1 * a_dim1] / b22;
00328 a22 = a[l1 + l1 * a_dim1] / b22;
00329 b12 = b[l + l1 * b_dim1] / b22;
00330 a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) / a21 +
00331 a12 - a11 * b12;
00332 a2 = a22 - a11 - a21 * b12 - (a33 - a11) - (a44 - a11) + a43 * b34;
00333 a3 = a[l1 + 1 + l1 * a_dim1] / b22;
00334 goto L160;
00335
00336 L155:
00337 a1 = 0.;
00338 a2 = 1.;
00339 a3 = 1.1605;
00340 L160:
00341 ++its;
00342 --itn;
00343 if (! (*matz)) {
00344 lor1 = ld;
00345 }
00346
00347 i__1 = na;
00348 for (k = l; k <= i__1; ++k) {
00349 notlas = k != na && ish == 2;
00350 k1 = k + 1;
00351 k2 = k + 2;
00352
00353 i__2 = k - 1;
00354 km1 = max(i__2,l);
00355
00356 i__2 = en, i__3 = k1 + ish;
00357 ll = min(i__2,i__3);
00358 if (notlas) {
00359 goto L190;
00360 }
00361
00362 if (k == l) {
00363 goto L170;
00364 }
00365 a1 = a[k + km1 * a_dim1];
00366 a2 = a[k1 + km1 * a_dim1];
00367 L170:
00368 s = abs(a1) + abs(a2);
00369 if (s == 0.) {
00370 goto L70;
00371 }
00372 u1 = a1 / s;
00373 u2 = a2 / s;
00374 d__1 = sqrt(u1 * u1 + u2 * u2);
00375 r__ = d_sign(&d__1, &u1);
00376 v1 = -(u1 + r__) / r__;
00377 v2 = -u2 / r__;
00378 u2 = v2 / v1;
00379
00380 i__2 = enorn;
00381 for (j = km1; j <= i__2; ++j) {
00382 t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1];
00383 a[k + j * a_dim1] += t * v1;
00384 a[k1 + j * a_dim1] += t * v2;
00385 t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1];
00386 b[k + j * b_dim1] += t * v1;
00387 b[k1 + j * b_dim1] += t * v2;
00388
00389 }
00390
00391 if (k != l) {
00392 a[k1 + km1 * a_dim1] = 0.;
00393 }
00394 goto L240;
00395
00396 L190:
00397 if (k == l) {
00398 goto L200;
00399 }
00400 a1 = a[k + km1 * a_dim1];
00401 a2 = a[k1 + km1 * a_dim1];
00402 a3 = a[k2 + km1 * a_dim1];
00403 L200:
00404 s = abs(a1) + abs(a2) + abs(a3);
00405 if (s == 0.) {
00406 goto L260;
00407 }
00408 u1 = a1 / s;
00409 u2 = a2 / s;
00410 u3 = a3 / s;
00411 d__1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
00412 r__ = d_sign(&d__1, &u1);
00413 v1 = -(u1 + r__) / r__;
00414 v2 = -u2 / r__;
00415 v3 = -u3 / r__;
00416 u2 = v2 / v1;
00417 u3 = v3 / v1;
00418
00419 i__2 = enorn;
00420 for (j = km1; j <= i__2; ++j) {
00421 t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1] + u3 * a[k2 + j *
00422 a_dim1];
00423 a[k + j * a_dim1] += t * v1;
00424 a[k1 + j * a_dim1] += t * v2;
00425 a[k2 + j * a_dim1] += t * v3;
00426 t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1] + u3 * b[k2 + j *
00427 b_dim1];
00428 b[k + j * b_dim1] += t * v1;
00429 b[k1 + j * b_dim1] += t * v2;
00430 b[k2 + j * b_dim1] += t * v3;
00431
00432 }
00433
00434 if (k == l) {
00435 goto L220;
00436 }
00437 a[k1 + km1 * a_dim1] = 0.;
00438 a[k2 + km1 * a_dim1] = 0.;
00439
00440 L220:
00441 s = (d__1 = b[k2 + k2 * b_dim1], abs(d__1)) + (d__2 = b[k2 + k1 *
00442 b_dim1], abs(d__2)) + (d__3 = b[k2 + k * b_dim1], abs(d__3));
00443 if (s == 0.) {
00444 goto L240;
00445 }
00446 u1 = b[k2 + k2 * b_dim1] / s;
00447 u2 = b[k2 + k1 * b_dim1] / s;
00448 u3 = b[k2 + k * b_dim1] / s;
00449 d__1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
00450 r__ = d_sign(&d__1, &u1);
00451 v1 = -(u1 + r__) / r__;
00452 v2 = -u2 / r__;
00453 v3 = -u3 / r__;
00454 u2 = v2 / v1;
00455 u3 = v3 / v1;
00456
00457 i__2 = ll;
00458 for (i__ = lor1; i__ <= i__2; ++i__) {
00459 t = a[i__ + k2 * a_dim1] + u2 * a[i__ + k1 * a_dim1] + u3 * a[i__
00460 + k * a_dim1];
00461 a[i__ + k2 * a_dim1] += t * v1;
00462 a[i__ + k1 * a_dim1] += t * v2;
00463 a[i__ + k * a_dim1] += t * v3;
00464 t = b[i__ + k2 * b_dim1] + u2 * b[i__ + k1 * b_dim1] + u3 * b[i__
00465 + k * b_dim1];
00466 b[i__ + k2 * b_dim1] += t * v1;
00467 b[i__ + k1 * b_dim1] += t * v2;
00468 b[i__ + k * b_dim1] += t * v3;
00469
00470 }
00471
00472 b[k2 + k * b_dim1] = 0.;
00473 b[k2 + k1 * b_dim1] = 0.;
00474 if (! (*matz)) {
00475 goto L240;
00476 }
00477
00478 i__2 = *n;
00479 for (i__ = 1; i__ <= i__2; ++i__) {
00480 t = z__[i__ + k2 * z_dim1] + u2 * z__[i__ + k1 * z_dim1] + u3 *
00481 z__[i__ + k * z_dim1];
00482 z__[i__ + k2 * z_dim1] += t * v1;
00483 z__[i__ + k1 * z_dim1] += t * v2;
00484 z__[i__ + k * z_dim1] += t * v3;
00485
00486 }
00487
00488 L240:
00489 s = (d__1 = b[k1 + k1 * b_dim1], abs(d__1)) + (d__2 = b[k1 + k *
00490 b_dim1], abs(d__2));
00491 if (s == 0.) {
00492 goto L260;
00493 }
00494 u1 = b[k1 + k1 * b_dim1] / s;
00495 u2 = b[k1 + k * b_dim1] / s;
00496 d__1 = sqrt(u1 * u1 + u2 * u2);
00497 r__ = d_sign(&d__1, &u1);
00498 v1 = -(u1 + r__) / r__;
00499 v2 = -u2 / r__;
00500 u2 = v2 / v1;
00501
00502 i__2 = ll;
00503 for (i__ = lor1; i__ <= i__2; ++i__) {
00504 t = a[i__ + k1 * a_dim1] + u2 * a[i__ + k * a_dim1];
00505 a[i__ + k1 * a_dim1] += t * v1;
00506 a[i__ + k * a_dim1] += t * v2;
00507 t = b[i__ + k1 * b_dim1] + u2 * b[i__ + k * b_dim1];
00508 b[i__ + k1 * b_dim1] += t * v1;
00509 b[i__ + k * b_dim1] += t * v2;
00510
00511 }
00512
00513 b[k1 + k * b_dim1] = 0.;
00514 if (! (*matz)) {
00515 goto L260;
00516 }
00517
00518 i__2 = *n;
00519 for (i__ = 1; i__ <= i__2; ++i__) {
00520 t = z__[i__ + k1 * z_dim1] + u2 * z__[i__ + k * z_dim1];
00521 z__[i__ + k1 * z_dim1] += t * v1;
00522 z__[i__ + k * z_dim1] += t * v2;
00523
00524 }
00525
00526 L260:
00527 ;
00528 }
00529
00530 goto L70;
00531
00532
00533 L1000:
00534 *ierr = en;
00535
00536 L1001:
00537 if (*n > 1) {
00538 b[*n + b_dim1] = epsb;
00539 }
00540 return 0;
00541 }
00542