Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_comlr.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int comlr_(integer *nm, integer *n, integer *low, integer *
00009 igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi,
00010 integer *ierr)
00011 {
00012
00013 integer hr_dim1, hr_offset, hi_dim1, hi_offset, i__1, i__2;
00014 doublereal d__1, d__2, d__3, d__4;
00015
00016
00017 extern int cdiv_(doublereal *, doublereal *, doublereal *
00018 , doublereal *, doublereal *, doublereal *);
00019 static integer i__, j, l, m, en, ll, mm;
00020 static doublereal si, ti, xi, yi, sr, tr, xr, yr;
00021 static integer im1;
00022 extern int csroot_(doublereal *, doublereal *,
00023 doublereal *, doublereal *);
00024 static integer mp1, itn, its;
00025 static doublereal zzi, zzr;
00026 static integer enm1;
00027 static doublereal tst1, tst2;
00028
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 --wi;
00087 --wr;
00088 hi_dim1 = *nm;
00089 hi_offset = hi_dim1 + 1;
00090 hi -= hi_offset;
00091 hr_dim1 = *nm;
00092 hr_offset = hr_dim1 + 1;
00093 hr -= hr_offset;
00094
00095
00096 *ierr = 0;
00097
00098 i__1 = *n;
00099 for (i__ = 1; i__ <= i__1; ++i__) {
00100 if (i__ >= *low && i__ <= *igh) {
00101 goto L200;
00102 }
00103 wr[i__] = hr[i__ + i__ * hr_dim1];
00104 wi[i__] = hi[i__ + i__ * hi_dim1];
00105 L200:
00106 ;
00107 }
00108
00109 en = *igh;
00110 tr = 0.;
00111 ti = 0.;
00112 itn = *n * 30;
00113
00114 L220:
00115 if (en < *low) {
00116 goto L1001;
00117 }
00118 its = 0;
00119 enm1 = en - 1;
00120
00121
00122 L240:
00123 i__1 = en;
00124 for (ll = *low; ll <= i__1; ++ll) {
00125 l = en + *low - ll;
00126 if (l == *low) {
00127 goto L300;
00128 }
00129 tst1 = (d__1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[
00130 l - 1 + (l - 1) * hi_dim1], abs(d__2)) + (d__3 = hr[l + l *
00131 hr_dim1], abs(d__3)) + (d__4 = hi[l + l * hi_dim1], abs(d__4))
00132 ;
00133 tst2 = tst1 + (d__1 = hr[l + (l - 1) * hr_dim1], abs(d__1)) + (d__2 =
00134 hi[l + (l - 1) * hi_dim1], abs(d__2));
00135 if (tst2 == tst1) {
00136 goto L300;
00137 }
00138
00139 }
00140
00141 L300:
00142 if (l == en) {
00143 goto L660;
00144 }
00145 if (itn == 0) {
00146 goto L1000;
00147 }
00148 if (its == 10 || its == 20) {
00149 goto L320;
00150 }
00151 sr = hr[en + en * hr_dim1];
00152 si = hi[en + en * hi_dim1];
00153 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en *
00154 hi_dim1] * hi[en + enm1 * hi_dim1];
00155 xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en *
00156 hi_dim1] * hr[en + enm1 * hr_dim1];
00157 if (xr == 0. && xi == 0.) {
00158 goto L340;
00159 }
00160 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
00161 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
00162
00163 d__2 = yr;
00164
00165 d__3 = yi;
00166 d__1 = d__2 * d__2 - d__3 * d__3 + xr;
00167 d__4 = yr * 2. * yi + xi;
00168 csroot_(&d__1, &d__4, &zzr, &zzi);
00169 if (yr * zzr + yi * zzi >= 0.) {
00170 goto L310;
00171 }
00172 zzr = -zzr;
00173 zzi = -zzi;
00174 L310:
00175 d__1 = yr + zzr;
00176 d__2 = yi + zzi;
00177 cdiv_(&xr, &xi, &d__1, &d__2, &xr, &xi);
00178 sr -= xr;
00179 si -= xi;
00180 goto L340;
00181
00182 L320:
00183 sr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hr[enm1 + (en
00184 - 2) * hr_dim1], abs(d__2));
00185 si = (d__1 = hi[en + enm1 * hi_dim1], abs(d__1)) + (d__2 = hi[enm1 + (en
00186 - 2) * hi_dim1], abs(d__2));
00187
00188 L340:
00189 i__1 = en;
00190 for (i__ = *low; i__ <= i__1; ++i__) {
00191 hr[i__ + i__ * hr_dim1] -= sr;
00192 hi[i__ + i__ * hi_dim1] -= si;
00193
00194 }
00195
00196 tr += sr;
00197 ti += si;
00198 ++its;
00199 --itn;
00200
00201
00202 xr = (d__1 = hr[enm1 + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[enm1 +
00203 enm1 * hi_dim1], abs(d__2));
00204 yr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hi[en + enm1 *
00205 hi_dim1], abs(d__2));
00206 zzr = (d__1 = hr[en + en * hr_dim1], abs(d__1)) + (d__2 = hi[en + en *
00207 hi_dim1], abs(d__2));
00208
00209 i__1 = enm1;
00210 for (mm = l; mm <= i__1; ++mm) {
00211 m = enm1 + l - mm;
00212 if (m == l) {
00213 goto L420;
00214 }
00215 yi = yr;
00216 yr = (d__1 = hr[m + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m + (
00217 m - 1) * hi_dim1], abs(d__2));
00218 xi = zzr;
00219 zzr = xr;
00220 xr = (d__1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[m
00221 - 1 + (m - 1) * hi_dim1], abs(d__2));
00222 tst1 = zzr / yi * (zzr + xr + xi);
00223 tst2 = tst1 + yr;
00224 if (tst2 == tst1) {
00225 goto L420;
00226 }
00227
00228 }
00229
00230 L420:
00231 mp1 = m + 1;
00232
00233 i__1 = en;
00234 for (i__ = mp1; i__ <= i__1; ++i__) {
00235 im1 = i__ - 1;
00236 xr = hr[im1 + im1 * hr_dim1];
00237 xi = hi[im1 + im1 * hi_dim1];
00238 yr = hr[i__ + im1 * hr_dim1];
00239 yi = hi[i__ + im1 * hi_dim1];
00240 if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
00241 goto L460;
00242 }
00243
00244 i__2 = en;
00245 for (j = im1; j <= i__2; ++j) {
00246 zzr = hr[im1 + j * hr_dim1];
00247 hr[im1 + j * hr_dim1] = hr[i__ + j * hr_dim1];
00248 hr[i__ + j * hr_dim1] = zzr;
00249 zzi = hi[im1 + j * hi_dim1];
00250 hi[im1 + j * hi_dim1] = hi[i__ + j * hi_dim1];
00251 hi[i__ + j * hi_dim1] = zzi;
00252
00253 }
00254
00255 cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
00256 wr[i__] = 1.;
00257 goto L480;
00258 L460:
00259 cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
00260 wr[i__] = -1.;
00261 L480:
00262 hr[i__ + im1 * hr_dim1] = zzr;
00263 hi[i__ + im1 * hi_dim1] = zzi;
00264
00265 i__2 = en;
00266 for (j = i__; j <= i__2; ++j) {
00267 hr[i__ + j * hr_dim1] = hr[i__ + j * hr_dim1] - zzr * hr[im1 + j *
00268 hr_dim1] + zzi * hi[im1 + j * hi_dim1];
00269 hi[i__ + j * hi_dim1] = hi[i__ + j * hi_dim1] - zzr * hi[im1 + j *
00270 hi_dim1] - zzi * hr[im1 + j * hr_dim1];
00271
00272 }
00273
00274
00275 }
00276
00277 i__1 = en;
00278 for (j = mp1; j <= i__1; ++j) {
00279 xr = hr[j + (j - 1) * hr_dim1];
00280 xi = hi[j + (j - 1) * hi_dim1];
00281 hr[j + (j - 1) * hr_dim1] = 0.;
00282 hi[j + (j - 1) * hi_dim1] = 0.;
00283
00284
00285 if (wr[j] <= 0.) {
00286 goto L580;
00287 }
00288
00289 i__2 = j;
00290 for (i__ = l; i__ <= i__2; ++i__) {
00291 zzr = hr[i__ + (j - 1) * hr_dim1];
00292 hr[i__ + (j - 1) * hr_dim1] = hr[i__ + j * hr_dim1];
00293 hr[i__ + j * hr_dim1] = zzr;
00294 zzi = hi[i__ + (j - 1) * hi_dim1];
00295 hi[i__ + (j - 1) * hi_dim1] = hi[i__ + j * hi_dim1];
00296 hi[i__ + j * hi_dim1] = zzi;
00297
00298 }
00299
00300 L580:
00301 i__2 = j;
00302 for (i__ = l; i__ <= i__2; ++i__) {
00303 hr[i__ + (j - 1) * hr_dim1] = hr[i__ + (j - 1) * hr_dim1] + xr *
00304 hr[i__ + j * hr_dim1] - xi * hi[i__ + j * hi_dim1];
00305 hi[i__ + (j - 1) * hi_dim1] = hi[i__ + (j - 1) * hi_dim1] + xr *
00306 hi[i__ + j * hi_dim1] + xi * hr[i__ + j * hr_dim1];
00307
00308 }
00309
00310
00311 }
00312
00313 goto L240;
00314
00315 L660:
00316 wr[en] = hr[en + en * hr_dim1] + tr;
00317 wi[en] = hi[en + en * hi_dim1] + ti;
00318 en = enm1;
00319 goto L220;
00320
00321
00322 L1000:
00323 *ierr = en;
00324 L1001:
00325 return 0;
00326 }
00327