00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int comlr2_(integer *nm, integer *n, integer *low, integer *
00009 igh, integer *int__, doublereal *hr, doublereal *hi, doublereal *wr,
00010 doublereal *wi, doublereal *zr, doublereal *zi, integer *ierr)
00011 {
00012
00013 integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset,
00014 zi_dim1, zi_offset, i__1, i__2, i__3;
00015 doublereal d__1, d__2, d__3, d__4;
00016
00017
00018 static integer iend;
00019 extern int cdiv_(doublereal *, doublereal *, doublereal *
00020 , doublereal *, doublereal *, doublereal *);
00021 static doublereal norm;
00022 static integer i__, j, k, l, m, ii, en, jj, ll, mm, nn;
00023 static doublereal si, ti, xi, yi, sr, tr, xr, yr;
00024 static integer im1;
00025 extern int csroot_(doublereal *, doublereal *,
00026 doublereal *, doublereal *);
00027 static integer ip1, mp1, itn, its;
00028 static doublereal zzi, zzr;
00029 static integer enm1;
00030 static doublereal tst1, tst2;
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
00100
00101
00102
00103
00104 zi_dim1 = *nm;
00105 zi_offset = zi_dim1 + 1;
00106 zi -= zi_offset;
00107 zr_dim1 = *nm;
00108 zr_offset = zr_dim1 + 1;
00109 zr -= zr_offset;
00110 --wi;
00111 --wr;
00112 hi_dim1 = *nm;
00113 hi_offset = hi_dim1 + 1;
00114 hi -= hi_offset;
00115 hr_dim1 = *nm;
00116 hr_offset = hr_dim1 + 1;
00117 hr -= hr_offset;
00118 --int__;
00119
00120
00121 *ierr = 0;
00122
00123 i__1 = *n;
00124 for (i__ = 1; i__ <= i__1; ++i__) {
00125
00126 i__2 = *n;
00127 for (j = 1; j <= i__2; ++j) {
00128 zr[i__ + j * zr_dim1] = 0.;
00129 zi[i__ + j * zi_dim1] = 0.;
00130 if (i__ == j) {
00131 zr[i__ + j * zr_dim1] = 1.;
00132 }
00133
00134 }
00135 }
00136
00137
00138 iend = *igh - *low - 1;
00139 if (iend <= 0) {
00140 goto L180;
00141 }
00142
00143 i__2 = iend;
00144 for (ii = 1; ii <= i__2; ++ii) {
00145 i__ = *igh - ii;
00146 ip1 = i__ + 1;
00147
00148 i__1 = *igh;
00149 for (k = ip1; k <= i__1; ++k) {
00150 zr[k + i__ * zr_dim1] = hr[k + (i__ - 1) * hr_dim1];
00151 zi[k + i__ * zi_dim1] = hi[k + (i__ - 1) * hi_dim1];
00152
00153 }
00154
00155 j = int__[i__];
00156 if (i__ == j) {
00157 goto L160;
00158 }
00159
00160 i__1 = *igh;
00161 for (k = i__; k <= i__1; ++k) {
00162 zr[i__ + k * zr_dim1] = zr[j + k * zr_dim1];
00163 zi[i__ + k * zi_dim1] = zi[j + k * zi_dim1];
00164 zr[j + k * zr_dim1] = 0.;
00165 zi[j + k * zi_dim1] = 0.;
00166
00167 }
00168
00169 zr[j + i__ * zr_dim1] = 1.;
00170 L160:
00171 ;
00172 }
00173
00174 L180:
00175 i__2 = *n;
00176 for (i__ = 1; i__ <= i__2; ++i__) {
00177 if (i__ >= *low && i__ <= *igh) {
00178 goto L200;
00179 }
00180 wr[i__] = hr[i__ + i__ * hr_dim1];
00181 wi[i__] = hi[i__ + i__ * hi_dim1];
00182 L200:
00183 ;
00184 }
00185
00186 en = *igh;
00187 tr = 0.;
00188 ti = 0.;
00189 itn = *n * 30;
00190
00191 L220:
00192 if (en < *low) {
00193 goto L680;
00194 }
00195 its = 0;
00196 enm1 = en - 1;
00197
00198
00199 L240:
00200 i__2 = en;
00201 for (ll = *low; ll <= i__2; ++ll) {
00202 l = en + *low - ll;
00203 if (l == *low) {
00204 goto L300;
00205 }
00206 tst1 = (d__1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[
00207 l - 1 + (l - 1) * hi_dim1], abs(d__2)) + (d__3 = hr[l + l *
00208 hr_dim1], abs(d__3)) + (d__4 = hi[l + l * hi_dim1], abs(d__4))
00209 ;
00210 tst2 = tst1 + (d__1 = hr[l + (l - 1) * hr_dim1], abs(d__1)) + (d__2 =
00211 hi[l + (l - 1) * hi_dim1], abs(d__2));
00212 if (tst2 == tst1) {
00213 goto L300;
00214 }
00215
00216 }
00217
00218 L300:
00219 if (l == en) {
00220 goto L660;
00221 }
00222 if (itn == 0) {
00223 goto L1000;
00224 }
00225 if (its == 10 || its == 20) {
00226 goto L320;
00227 }
00228 sr = hr[en + en * hr_dim1];
00229 si = hi[en + en * hi_dim1];
00230 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en *
00231 hi_dim1] * hi[en + enm1 * hi_dim1];
00232 xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en *
00233 hi_dim1] * hr[en + enm1 * hr_dim1];
00234 if (xr == 0. && xi == 0.) {
00235 goto L340;
00236 }
00237 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
00238 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
00239
00240 d__2 = yr;
00241
00242 d__3 = yi;
00243 d__1 = d__2 * d__2 - d__3 * d__3 + xr;
00244 d__4 = yr * 2. * yi + xi;
00245 csroot_(&d__1, &d__4, &zzr, &zzi);
00246 if (yr * zzr + yi * zzi >= 0.) {
00247 goto L310;
00248 }
00249 zzr = -zzr;
00250 zzi = -zzi;
00251 L310:
00252 d__1 = yr + zzr;
00253 d__2 = yi + zzi;
00254 cdiv_(&xr, &xi, &d__1, &d__2, &xr, &xi);
00255 sr -= xr;
00256 si -= xi;
00257 goto L340;
00258
00259 L320:
00260 sr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hr[enm1 + (en
00261 - 2) * hr_dim1], abs(d__2));
00262 si = (d__1 = hi[en + enm1 * hi_dim1], abs(d__1)) + (d__2 = hi[enm1 + (en
00263 - 2) * hi_dim1], abs(d__2));
00264
00265 L340:
00266 i__2 = en;
00267 for (i__ = *low; i__ <= i__2; ++i__) {
00268 hr[i__ + i__ * hr_dim1] -= sr;
00269 hi[i__ + i__ * hi_dim1] -= si;
00270
00271 }
00272
00273 tr += sr;
00274 ti += si;
00275 ++its;
00276 --itn;
00277
00278
00279 xr = (d__1 = hr[enm1 + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[enm1 +
00280 enm1 * hi_dim1], abs(d__2));
00281 yr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[en + enm1 *
00282 hi_dim1], abs(d__2));
00283 zzr = (d__1 = hr[en + en * hr_dim1], abs(d__1)) + (d__2 = hi[en + en *
00284 hi_dim1], abs(d__2));
00285
00286 i__2 = enm1;
00287 for (mm = l; mm <= i__2; ++mm) {
00288 m = enm1 + l - mm;
00289 if (m == l) {
00290 goto L420;
00291 }
00292 yi = yr;
00293 yr = (d__1 = hr[m + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m + (
00294 m - 1) * hi_dim1], abs(d__2));
00295 xi = zzr;
00296 zzr = xr;
00297 xr = (d__1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m
00298 - 1 + (m - 1) * hi_dim1], abs(d__2));
00299 tst1 = zzr / yi * (zzr + xr + xi);
00300 tst2 = tst1 + yr;
00301 if (tst2 == tst1) {
00302 goto L420;
00303 }
00304
00305 }
00306
00307 L420:
00308 mp1 = m + 1;
00309
00310 i__2 = en;
00311 for (i__ = mp1; i__ <= i__2; ++i__) {
00312 im1 = i__ - 1;
00313 xr = hr[im1 + im1 * hr_dim1];
00314 xi = hi[im1 + im1 * hi_dim1];
00315 yr = hr[i__ + im1 * hr_dim1];
00316 yi = hi[i__ + im1 * hi_dim1];
00317 if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
00318 goto L460;
00319 }
00320
00321 i__1 = *n;
00322 for (j = im1; j <= i__1; ++j) {
00323 zzr = hr[im1 + j * hr_dim1];
00324 hr[im1 + j * hr_dim1] = hr[i__ + j * hr_dim1];
00325 hr[i__ + j * hr_dim1] = zzr;
00326 zzi = hi[im1 + j * hi_dim1];
00327 hi[im1 + j * hi_dim1] = hi[i__ + j * hi_dim1];
00328 hi[i__ + j * hi_dim1] = zzi;
00329
00330 }
00331
00332 cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
00333 wr[i__] = 1.;
00334 goto L480;
00335 L460:
00336 cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
00337 wr[i__] = -1.;
00338 L480:
00339 hr[i__ + im1 * hr_dim1] = zzr;
00340 hi[i__ + im1 * hi_dim1] = zzi;
00341
00342 i__1 = *n;
00343 for (j = i__; j <= i__1; ++j) {
00344 hr[i__ + j * hr_dim1] = hr[i__ + j * hr_dim1] - zzr * hr[im1 + j *
00345 hr_dim1] + zzi * hi[im1 + j * hi_dim1];
00346 hi[i__ + j * hi_dim1] = hi[i__ + j * hi_dim1] - zzr * hi[im1 + j *
00347 hi_dim1] - zzi * hr[im1 + j * hr_dim1];
00348
00349 }
00350
00351
00352 }
00353
00354 i__2 = en;
00355 for (j = mp1; j <= i__2; ++j) {
00356 xr = hr[j + (j - 1) * hr_dim1];
00357 xi = hi[j + (j - 1) * hi_dim1];
00358 hr[j + (j - 1) * hr_dim1] = 0.;
00359 hi[j + (j - 1) * hi_dim1] = 0.;
00360
00361
00362 if (wr[j] <= 0.) {
00363 goto L580;
00364 }
00365
00366 i__1 = j;
00367 for (i__ = 1; i__ <= i__1; ++i__) {
00368 zzr = hr[i__ + (j - 1) * hr_dim1];
00369 hr[i__ + (j - 1) * hr_dim1] = hr[i__ + j * hr_dim1];
00370 hr[i__ + j * hr_dim1] = zzr;
00371 zzi = hi[i__ + (j - 1) * hi_dim1];
00372 hi[i__ + (j - 1) * hi_dim1] = hi[i__ + j * hi_dim1];
00373 hi[i__ + j * hi_dim1] = zzi;
00374
00375 }
00376
00377 i__1 = *igh;
00378 for (i__ = *low; i__ <= i__1; ++i__) {
00379 zzr = zr[i__ + (j - 1) * zr_dim1];
00380 zr[i__ + (j - 1) * zr_dim1] = zr[i__ + j * zr_dim1];
00381 zr[i__ + j * zr_dim1] = zzr;
00382 zzi = zi[i__ + (j - 1) * zi_dim1];
00383 zi[i__ + (j - 1) * zi_dim1] = zi[i__ + j * zi_dim1];
00384 zi[i__ + j * zi_dim1] = zzi;
00385
00386 }
00387
00388 L580:
00389 i__1 = j;
00390 for (i__ = 1; i__ <= i__1; ++i__) {
00391 hr[i__ + (j - 1) * hr_dim1] = hr[i__ + (j - 1) * hr_dim1] + xr *
00392 hr[i__ + j * hr_dim1] - xi * hi[i__ + j * hi_dim1];
00393 hi[i__ + (j - 1) * hi_dim1] = hi[i__ + (j - 1) * hi_dim1] + xr *
00394 hi[i__ + j * hi_dim1] + xi * hr[i__ + j * hr_dim1];
00395
00396 }
00397
00398 i__1 = *igh;
00399 for (i__ = *low; i__ <= i__1; ++i__) {
00400 zr[i__ + (j - 1) * zr_dim1] = zr[i__ + (j - 1) * zr_dim1] + xr *
00401 zr[i__ + j * zr_dim1] - xi * zi[i__ + j * zi_dim1];
00402 zi[i__ + (j - 1) * zi_dim1] = zi[i__ + (j - 1) * zi_dim1] + xr *
00403 zi[i__ + j * zi_dim1] + xi * zr[i__ + j * zr_dim1];
00404
00405 }
00406
00407
00408 }
00409
00410 goto L240;
00411
00412 L660:
00413 hr[en + en * hr_dim1] += tr;
00414 wr[en] = hr[en + en * hr_dim1];
00415 hi[en + en * hi_dim1] += ti;
00416 wi[en] = hi[en + en * hi_dim1];
00417 en = enm1;
00418 goto L220;
00419
00420
00421 L680:
00422 norm = 0.;
00423
00424 i__2 = *n;
00425 for (i__ = 1; i__ <= i__2; ++i__) {
00426
00427 i__1 = *n;
00428 for (j = i__; j <= i__1; ++j) {
00429 tr = (d__1 = hr[i__ + j * hr_dim1], abs(d__1)) + (d__2 = hi[i__ +
00430 j * hi_dim1], abs(d__2));
00431 if (tr > norm) {
00432 norm = tr;
00433 }
00434
00435 }
00436 }
00437
00438 hr[hr_dim1 + 1] = norm;
00439 if (*n == 1 || norm == 0.) {
00440 goto L1001;
00441 }
00442
00443 i__1 = *n;
00444 for (nn = 2; nn <= i__1; ++nn) {
00445 en = *n + 2 - nn;
00446 xr = wr[en];
00447 xi = wi[en];
00448 hr[en + en * hr_dim1] = 1.;
00449 hi[en + en * hi_dim1] = 0.;
00450 enm1 = en - 1;
00451
00452 i__2 = enm1;
00453 for (ii = 1; ii <= i__2; ++ii) {
00454 i__ = en - ii;
00455 zzr = 0.;
00456 zzi = 0.;
00457 ip1 = i__ + 1;
00458
00459 i__3 = en;
00460 for (j = ip1; j <= i__3; ++j) {
00461 zzr = zzr + hr[i__ + j * hr_dim1] * hr[j + en * hr_dim1] - hi[
00462 i__ + j * hi_dim1] * hi[j + en * hi_dim1];
00463 zzi = zzi + hr[i__ + j * hr_dim1] * hi[j + en * hi_dim1] + hi[
00464 i__ + j * hi_dim1] * hr[j + en * hr_dim1];
00465
00466 }
00467
00468 yr = xr - wr[i__];
00469 yi = xi - wi[i__];
00470 if (yr != 0. || yi != 0.) {
00471 goto L765;
00472 }
00473 tst1 = norm;
00474 yr = tst1;
00475 L760:
00476 yr *= .01;
00477 tst2 = norm + yr;
00478 if (tst2 > tst1) {
00479 goto L760;
00480 }
00481 L765:
00482 cdiv_(&zzr, &zzi, &yr, &yi, &hr[i__ + en * hr_dim1], &hi[i__ + en
00483 * hi_dim1]);
00484
00485 tr = (d__1 = hr[i__ + en * hr_dim1], abs(d__1)) + (d__2 = hi[i__
00486 + en * hi_dim1], abs(d__2));
00487 if (tr == 0.) {
00488 goto L780;
00489 }
00490 tst1 = tr;
00491 tst2 = tst1 + 1. / tst1;
00492 if (tst2 > tst1) {
00493 goto L780;
00494 }
00495 i__3 = en;
00496 for (j = i__; j <= i__3; ++j) {
00497 hr[j + en * hr_dim1] /= tr;
00498 hi[j + en * hi_dim1] /= tr;
00499
00500 }
00501
00502 L780:
00503 ;
00504 }
00505
00506
00507 }
00508
00509 enm1 = *n - 1;
00510
00511 i__1 = enm1;
00512 for (i__ = 1; i__ <= i__1; ++i__) {
00513 if (i__ >= *low && i__ <= *igh) {
00514 goto L840;
00515 }
00516 ip1 = i__ + 1;
00517
00518 i__2 = *n;
00519 for (j = ip1; j <= i__2; ++j) {
00520 zr[i__ + j * zr_dim1] = hr[i__ + j * hr_dim1];
00521 zi[i__ + j * zi_dim1] = hi[i__ + j * hi_dim1];
00522
00523 }
00524
00525 L840:
00526 ;
00527 }
00528
00529
00530
00531 i__1 = enm1;
00532 for (jj = *low; jj <= i__1; ++jj) {
00533 j = *n + *low - jj;
00534 m = min(j,*igh);
00535
00536 i__2 = *igh;
00537 for (i__ = *low; i__ <= i__2; ++i__) {
00538 zzr = 0.;
00539 zzi = 0.;
00540
00541 i__3 = m;
00542 for (k = *low; k <= i__3; ++k) {
00543 zzr = zzr + zr[i__ + k * zr_dim1] * hr[k + j * hr_dim1] - zi[
00544 i__ + k * zi_dim1] * hi[k + j * hi_dim1];
00545 zzi = zzi + zr[i__ + k * zr_dim1] * hi[k + j * hi_dim1] + zi[
00546 i__ + k * zi_dim1] * hr[k + j * hr_dim1];
00547
00548 }
00549
00550 zr[i__ + j * zr_dim1] = zzr;
00551 zi[i__ + j * zi_dim1] = zzi;
00552
00553 }
00554 }
00555
00556 goto L1001;
00557
00558
00559 L1000:
00560 *ierr = en;
00561 L1001:
00562 return 0;
00563 }
00564