Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_invit.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int invit_(integer *nm, integer *n, doublereal *a,
00009 doublereal *wr, doublereal *wi, logical *select, integer *mm, integer
00010 *m, doublereal *z__, integer *ierr, doublereal *rm1, doublereal *rv1,
00011 doublereal *rv2)
00012 {
00013
00014 integer a_dim1, a_offset, z_dim1, z_offset, rm1_dim1, rm1_offset, i__1,
00015 i__2, i__3;
00016 doublereal d__1, d__2;
00017
00018
00019 double sqrt(doublereal);
00020
00021
00022 extern int cdiv_(doublereal *, doublereal *, doublereal *
00023 , doublereal *, doublereal *, doublereal *);
00024 static doublereal norm;
00025 static integer i__, j, k, l, s;
00026 static doublereal t, w, x, y;
00027 static integer n1;
00028 static doublereal normv;
00029 static integer ii;
00030 static doublereal ilambd;
00031 static integer ip, mp, ns, uk;
00032 static doublereal rlambd;
00033 extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal
00034 *);
00035 static integer km1, ip1;
00036 static doublereal growto, ukroot;
00037 static integer its;
00038 static doublereal eps3;
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
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127 --rv2;
00128 --rv1;
00129 rm1_dim1 = *n;
00130 rm1_offset = rm1_dim1 + 1;
00131 rm1 -= rm1_offset;
00132 --select;
00133 --wi;
00134 --wr;
00135 a_dim1 = *nm;
00136 a_offset = a_dim1 + 1;
00137 a -= a_offset;
00138 z_dim1 = *nm;
00139 z_offset = z_dim1 + 1;
00140 z__ -= z_offset;
00141
00142
00143 *ierr = 0;
00144 uk = 0;
00145 s = 1;
00146
00147
00148
00149 ip = 0;
00150 n1 = *n - 1;
00151
00152 i__1 = *n;
00153 for (k = 1; k <= i__1; ++k) {
00154 if (wi[k] == 0. || ip < 0) {
00155 goto L100;
00156 }
00157 ip = 1;
00158 if (select[k] && select[k + 1]) {
00159 select[k + 1] = FALSE_;
00160 }
00161 L100:
00162 if (! select[k]) {
00163 goto L960;
00164 }
00165 if (wi[k] != 0.) {
00166 ++s;
00167 }
00168 if (s > *mm) {
00169 goto L1000;
00170 }
00171 if (uk >= k) {
00172 goto L200;
00173 }
00174
00175 i__2 = *n;
00176 for (uk = k; uk <= i__2; ++uk) {
00177 if (uk == *n) {
00178 goto L140;
00179 }
00180 if (a[uk + 1 + uk * a_dim1] == 0.) {
00181 goto L140;
00182 }
00183
00184 }
00185
00186
00187 L140:
00188 norm = 0.;
00189 mp = 1;
00190
00191 i__2 = uk;
00192 for (i__ = 1; i__ <= i__2; ++i__) {
00193 x = 0.;
00194
00195 i__3 = uk;
00196 for (j = mp; j <= i__3; ++j) {
00197
00198 x += (d__1 = a[i__ + j * a_dim1], abs(d__1));
00199 }
00200
00201 if (x > norm) {
00202 norm = x;
00203 }
00204 mp = i__;
00205
00206 }
00207
00208
00209 if (norm == 0.) {
00210 norm = 1.;
00211 }
00212 eps3 = epslon_(&norm);
00213
00214 ukroot = (doublereal) uk;
00215 ukroot = sqrt(ukroot);
00216 growto = .1 / ukroot;
00217 L200:
00218 rlambd = wr[k];
00219 ilambd = wi[k];
00220 if (k == 1) {
00221 goto L280;
00222 }
00223 km1 = k - 1;
00224 goto L240;
00225
00226
00227 L220:
00228 rlambd += eps3;
00229
00230 L240:
00231 i__2 = km1;
00232 for (ii = 1; ii <= i__2; ++ii) {
00233 i__ = k - ii;
00234 if (select[i__] && (d__1 = wr[i__] - rlambd, abs(d__1)) < eps3 &&
00235 (d__2 = wi[i__] - ilambd, abs(d__2)) < eps3) {
00236 goto L220;
00237 }
00238
00239 }
00240
00241 wr[k] = rlambd;
00242
00243 ip1 = k + ip;
00244 wr[ip1] = rlambd;
00245
00246
00247 L280:
00248 mp = 1;
00249
00250 i__2 = uk;
00251 for (i__ = 1; i__ <= i__2; ++i__) {
00252
00253 i__3 = uk;
00254 for (j = mp; j <= i__3; ++j) {
00255
00256 rm1[j + i__ * rm1_dim1] = a[i__ + j * a_dim1];
00257 }
00258
00259 rm1[i__ + i__ * rm1_dim1] -= rlambd;
00260 mp = i__;
00261 rv1[i__] = eps3;
00262
00263 }
00264
00265 its = 0;
00266 if (ilambd != 0.) {
00267 goto L520;
00268 }
00269
00270
00271
00272 if (uk == 1) {
00273 goto L420;
00274 }
00275
00276 i__2 = uk;
00277 for (i__ = 2; i__ <= i__2; ++i__) {
00278 mp = i__ - 1;
00279 if ((d__1 = rm1[mp + i__ * rm1_dim1], abs(d__1)) <= (d__2 = rm1[
00280 mp + mp * rm1_dim1], abs(d__2))) {
00281 goto L360;
00282 }
00283
00284 i__3 = uk;
00285 for (j = mp; j <= i__3; ++j) {
00286 y = rm1[j + i__ * rm1_dim1];
00287 rm1[j + i__ * rm1_dim1] = rm1[j + mp * rm1_dim1];
00288 rm1[j + mp * rm1_dim1] = y;
00289
00290 }
00291
00292 L360:
00293 if (rm1[mp + mp * rm1_dim1] == 0.) {
00294 rm1[mp + mp * rm1_dim1] = eps3;
00295 }
00296 x = rm1[mp + i__ * rm1_dim1] / rm1[mp + mp * rm1_dim1];
00297 if (x == 0.) {
00298 goto L400;
00299 }
00300
00301 i__3 = uk;
00302 for (j = i__; j <= i__3; ++j) {
00303
00304 rm1[j + i__ * rm1_dim1] -= x * rm1[j + mp * rm1_dim1];
00305 }
00306
00307 L400:
00308 ;
00309 }
00310
00311 L420:
00312 if (rm1[uk + uk * rm1_dim1] == 0.) {
00313 rm1[uk + uk * rm1_dim1] = eps3;
00314 }
00315
00316
00317 L440:
00318 i__2 = uk;
00319 for (ii = 1; ii <= i__2; ++ii) {
00320 i__ = uk + 1 - ii;
00321 y = rv1[i__];
00322 if (i__ == uk) {
00323 goto L480;
00324 }
00325 ip1 = i__ + 1;
00326
00327 i__3 = uk;
00328 for (j = ip1; j <= i__3; ++j) {
00329
00330 y -= rm1[j + i__ * rm1_dim1] * rv1[j];
00331 }
00332
00333 L480:
00334 rv1[i__] = y / rm1[i__ + i__ * rm1_dim1];
00335
00336 }
00337
00338 goto L740;
00339
00340
00341
00342
00343
00344 L520:
00345 ns = *n - s;
00346 z__[(s - 1) * z_dim1 + 1] = -ilambd;
00347 z__[s * z_dim1 + 1] = 0.;
00348 if (*n == 2) {
00349 goto L550;
00350 }
00351 rm1[rm1_dim1 * 3 + 1] = -ilambd;
00352 z__[(s - 1) * z_dim1 + 1] = 0.;
00353 if (*n == 3) {
00354 goto L550;
00355 }
00356
00357 i__2 = *n;
00358 for (i__ = 4; i__ <= i__2; ++i__) {
00359
00360 rm1[i__ * rm1_dim1 + 1] = 0.;
00361 }
00362
00363 L550:
00364 i__2 = uk;
00365 for (i__ = 2; i__ <= i__2; ++i__) {
00366 mp = i__ - 1;
00367 w = rm1[mp + i__ * rm1_dim1];
00368 if (i__ < *n) {
00369 t = rm1[mp + (i__ + 1) * rm1_dim1];
00370 }
00371 if (i__ == *n) {
00372 t = z__[mp + (s - 1) * z_dim1];
00373 }
00374 x = rm1[mp + mp * rm1_dim1] * rm1[mp + mp * rm1_dim1] + t * t;
00375 if (w * w <= x) {
00376 goto L580;
00377 }
00378 x = rm1[mp + mp * rm1_dim1] / w;
00379 y = t / w;
00380 rm1[mp + mp * rm1_dim1] = w;
00381 if (i__ < *n) {
00382 rm1[mp + (i__ + 1) * rm1_dim1] = 0.;
00383 }
00384 if (i__ == *n) {
00385 z__[mp + (s - 1) * z_dim1] = 0.;
00386 }
00387
00388 i__3 = uk;
00389 for (j = i__; j <= i__3; ++j) {
00390 w = rm1[j + i__ * rm1_dim1];
00391 rm1[j + i__ * rm1_dim1] = rm1[j + mp * rm1_dim1] - x * w;
00392 rm1[j + mp * rm1_dim1] = w;
00393 if (j < n1) {
00394 goto L555;
00395 }
00396 l = j - ns;
00397 z__[i__ + l * z_dim1] = z__[mp + l * z_dim1] - y * w;
00398 z__[mp + l * z_dim1] = 0.;
00399 goto L560;
00400 L555:
00401 rm1[i__ + (j + 2) * rm1_dim1] = rm1[mp + (j + 2) * rm1_dim1]
00402 - y * w;
00403 rm1[mp + (j + 2) * rm1_dim1] = 0.;
00404 L560:
00405 ;
00406 }
00407
00408 rm1[i__ + i__ * rm1_dim1] -= y * ilambd;
00409 if (i__ < n1) {
00410 goto L570;
00411 }
00412 l = i__ - ns;
00413 z__[mp + l * z_dim1] = -ilambd;
00414 z__[i__ + l * z_dim1] += x * ilambd;
00415 goto L640;
00416 L570:
00417 rm1[mp + (i__ + 2) * rm1_dim1] = -ilambd;
00418 rm1[i__ + (i__ + 2) * rm1_dim1] += x * ilambd;
00419 goto L640;
00420 L580:
00421 if (x != 0.) {
00422 goto L600;
00423 }
00424 rm1[mp + mp * rm1_dim1] = eps3;
00425 if (i__ < *n) {
00426 rm1[mp + (i__ + 1) * rm1_dim1] = 0.;
00427 }
00428 if (i__ == *n) {
00429 z__[mp + (s - 1) * z_dim1] = 0.;
00430 }
00431 t = 0.;
00432 x = eps3 * eps3;
00433 L600:
00434 w /= x;
00435 x = rm1[mp + mp * rm1_dim1] * w;
00436 y = -t * w;
00437
00438 i__3 = uk;
00439 for (j = i__; j <= i__3; ++j) {
00440 if (j < n1) {
00441 goto L610;
00442 }
00443 l = j - ns;
00444 t = z__[mp + l * z_dim1];
00445 z__[i__ + l * z_dim1] = -x * t - y * rm1[j + mp * rm1_dim1];
00446 goto L615;
00447 L610:
00448 t = rm1[mp + (j + 2) * rm1_dim1];
00449 rm1[i__ + (j + 2) * rm1_dim1] = -x * t - y * rm1[j + mp *
00450 rm1_dim1];
00451 L615:
00452 rm1[j + i__ * rm1_dim1] = rm1[j + i__ * rm1_dim1] - x * rm1[j
00453 + mp * rm1_dim1] + y * t;
00454
00455 }
00456
00457 if (i__ < n1) {
00458 goto L630;
00459 }
00460 l = i__ - ns;
00461 z__[i__ + l * z_dim1] -= ilambd;
00462 goto L640;
00463 L630:
00464 rm1[i__ + (i__ + 2) * rm1_dim1] -= ilambd;
00465 L640:
00466 ;
00467 }
00468
00469 if (uk < n1) {
00470 goto L650;
00471 }
00472 l = uk - ns;
00473 t = z__[uk + l * z_dim1];
00474 goto L655;
00475 L650:
00476 t = rm1[uk + (uk + 2) * rm1_dim1];
00477 L655:
00478 if (rm1[uk + uk * rm1_dim1] == 0. && t == 0.) {
00479 rm1[uk + uk * rm1_dim1] = eps3;
00480 }
00481
00482
00483 L660:
00484 i__2 = uk;
00485 for (ii = 1; ii <= i__2; ++ii) {
00486 i__ = uk + 1 - ii;
00487 x = rv1[i__];
00488 y = 0.;
00489 if (i__ == uk) {
00490 goto L700;
00491 }
00492 ip1 = i__ + 1;
00493
00494 i__3 = uk;
00495 for (j = ip1; j <= i__3; ++j) {
00496 if (j < n1) {
00497 goto L670;
00498 }
00499 l = j - ns;
00500 t = z__[i__ + l * z_dim1];
00501 goto L675;
00502 L670:
00503 t = rm1[i__ + (j + 2) * rm1_dim1];
00504 L675:
00505 x = x - rm1[j + i__ * rm1_dim1] * rv1[j] + t * rv2[j];
00506 y = y - rm1[j + i__ * rm1_dim1] * rv2[j] - t * rv1[j];
00507
00508 }
00509
00510 L700:
00511 if (i__ < n1) {
00512 goto L710;
00513 }
00514 l = i__ - ns;
00515 t = z__[i__ + l * z_dim1];
00516 goto L715;
00517 L710:
00518 t = rm1[i__ + (i__ + 2) * rm1_dim1];
00519 L715:
00520 cdiv_(&x, &y, &rm1[i__ + i__ * rm1_dim1], &t, &rv1[i__], &rv2[i__]
00521 );
00522
00523 }
00524
00525
00526 L740:
00527 ++its;
00528 norm = 0.;
00529 normv = 0.;
00530
00531 i__2 = uk;
00532 for (i__ = 1; i__ <= i__2; ++i__) {
00533 if (ilambd == 0.) {
00534 x = (d__1 = rv1[i__], abs(d__1));
00535 }
00536 if (ilambd != 0.) {
00537 x = pythag_(&rv1[i__], &rv2[i__]);
00538 }
00539 if (normv >= x) {
00540 goto L760;
00541 }
00542 normv = x;
00543 j = i__;
00544 L760:
00545 norm += x;
00546
00547 }
00548
00549 if (norm < growto) {
00550 goto L840;
00551 }
00552
00553 x = rv1[j];
00554 if (ilambd == 0.) {
00555 x = 1. / x;
00556 }
00557 if (ilambd != 0.) {
00558 y = rv2[j];
00559 }
00560
00561 i__2 = uk;
00562 for (i__ = 1; i__ <= i__2; ++i__) {
00563 if (ilambd != 0.) {
00564 goto L800;
00565 }
00566 z__[i__ + s * z_dim1] = rv1[i__] * x;
00567 goto L820;
00568 L800:
00569 cdiv_(&rv1[i__], &rv2[i__], &x, &y, &z__[i__ + (s - 1) * z_dim1],
00570 &z__[i__ + s * z_dim1]);
00571 L820:
00572 ;
00573 }
00574
00575 if (uk == *n) {
00576 goto L940;
00577 }
00578 j = uk + 1;
00579 goto L900;
00580
00581
00582 L840:
00583 if (its >= uk) {
00584 goto L880;
00585 }
00586 x = ukroot;
00587 y = eps3 / (x + 1.);
00588 rv1[1] = eps3;
00589
00590 i__2 = uk;
00591 for (i__ = 2; i__ <= i__2; ++i__) {
00592
00593 rv1[i__] = y;
00594 }
00595
00596 j = uk - its + 1;
00597 rv1[j] -= eps3 * x;
00598 if (ilambd == 0.) {
00599 goto L440;
00600 }
00601 goto L660;
00602
00603 L880:
00604 j = 1;
00605 *ierr = -k;
00606
00607
00608 L900:
00609 i__2 = *n;
00610 for (i__ = j; i__ <= i__2; ++i__) {
00611 z__[i__ + s * z_dim1] = 0.;
00612 if (ilambd != 0.) {
00613 z__[i__ + (s - 1) * z_dim1] = 0.;
00614 }
00615
00616 }
00617
00618 L940:
00619 ++s;
00620 L960:
00621 if (ip == -1) {
00622 ip = 0;
00623 }
00624 if (ip == 1) {
00625 ip = -1;
00626 }
00627
00628 }
00629
00630 goto L1001;
00631
00632
00633 L1000:
00634 if (*ierr != 0) {
00635 *ierr -= *n;
00636 }
00637 if (*ierr == 0) {
00638 *ierr = -((*n << 1) + 1);
00639 }
00640 L1001:
00641 *m = s - 1 - abs(ip);
00642 return 0;
00643 }
00644