Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_hqr.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int hqr_(integer *nm, integer *n, integer *low, integer *igh,
00009 doublereal *h__, doublereal *wr, doublereal *wi, integer *ierr)
00010 {
00011
00012 integer h_dim1, h_offset, i__1, i__2, i__3;
00013 doublereal d__1, d__2;
00014
00015
00016 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00017
00018
00019 static doublereal norm;
00020 static integer i__, j, k, l, m;
00021 static doublereal p, q, r__, s, t, w, x, y;
00022 static integer na, en, ll, mm;
00023 static doublereal zz;
00024 static logical notlas;
00025 static integer mp2, itn, its, enm2;
00026 static doublereal tst1, tst2;
00027
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 --wi;
00084 --wr;
00085 h_dim1 = *nm;
00086 h_offset = h_dim1 + 1;
00087 h__ -= h_offset;
00088
00089
00090 *ierr = 0;
00091 norm = 0.;
00092 k = 1;
00093
00094
00095 i__1 = *n;
00096 for (i__ = 1; i__ <= i__1; ++i__) {
00097
00098 i__2 = *n;
00099 for (j = k; j <= i__2; ++j) {
00100
00101 norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1));
00102 }
00103
00104 k = i__;
00105 if (i__ >= *low && i__ <= *igh) {
00106 goto L50;
00107 }
00108 wr[i__] = h__[i__ + i__ * h_dim1];
00109 wi[i__] = 0.;
00110 L50:
00111 ;
00112 }
00113
00114 en = *igh;
00115 t = 0.;
00116 itn = *n * 30;
00117
00118 L60:
00119 if (en < *low) {
00120 goto L1001;
00121 }
00122 its = 0;
00123 na = en - 1;
00124 enm2 = na - 1;
00125
00126
00127 L70:
00128 i__1 = en;
00129 for (ll = *low; ll <= i__1; ++ll) {
00130 l = en + *low - ll;
00131 if (l == *low) {
00132 goto L100;
00133 }
00134 s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l
00135 + l * h_dim1], abs(d__2));
00136 if (s == 0.) {
00137 s = norm;
00138 }
00139 tst1 = s;
00140 tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1));
00141 if (tst2 == tst1) {
00142 goto L100;
00143 }
00144
00145 }
00146
00147 L100:
00148 x = h__[en + en * h_dim1];
00149 if (l == en) {
00150 goto L270;
00151 }
00152 y = h__[na + na * h_dim1];
00153 w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
00154 if (l == na) {
00155 goto L280;
00156 }
00157 if (itn == 0) {
00158 goto L1000;
00159 }
00160 if (its != 10 && its != 20) {
00161 goto L130;
00162 }
00163
00164 t += x;
00165
00166 i__1 = en;
00167 for (i__ = *low; i__ <= i__1; ++i__) {
00168
00169 h__[i__ + i__ * h_dim1] -= x;
00170 }
00171
00172 s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 *
00173 h_dim1], abs(d__2));
00174 x = s * .75;
00175 y = x;
00176 w = s * -.4375 * s;
00177 L130:
00178 ++its;
00179 --itn;
00180
00181
00182
00183 i__1 = enm2;
00184 for (mm = l; mm <= i__1; ++mm) {
00185 m = enm2 + l - mm;
00186 zz = h__[m + m * h_dim1];
00187 r__ = x - zz;
00188 s = y - zz;
00189 p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) *
00190 h_dim1];
00191 q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
00192 r__ = h__[m + 2 + (m + 1) * h_dim1];
00193 s = abs(p) + abs(q) + abs(r__);
00194 p /= s;
00195 q /= s;
00196 r__ /= s;
00197 if (m == l) {
00198 goto L150;
00199 }
00200 tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) +
00201 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
00202 tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q)
00203 + abs(r__));
00204 if (tst2 == tst1) {
00205 goto L150;
00206 }
00207
00208 }
00209
00210 L150:
00211 mp2 = m + 2;
00212
00213 i__1 = en;
00214 for (i__ = mp2; i__ <= i__1; ++i__) {
00215 h__[i__ + (i__ - 2) * h_dim1] = 0.;
00216 if (i__ == mp2) {
00217 goto L160;
00218 }
00219 h__[i__ + (i__ - 3) * h_dim1] = 0.;
00220 L160:
00221 ;
00222 }
00223
00224
00225 i__1 = na;
00226 for (k = m; k <= i__1; ++k) {
00227 notlas = k != na;
00228 if (k == m) {
00229 goto L170;
00230 }
00231 p = h__[k + (k - 1) * h_dim1];
00232 q = h__[k + 1 + (k - 1) * h_dim1];
00233 r__ = 0.;
00234 if (notlas) {
00235 r__ = h__[k + 2 + (k - 1) * h_dim1];
00236 }
00237 x = abs(p) + abs(q) + abs(r__);
00238 if (x == 0.) {
00239 goto L260;
00240 }
00241 p /= x;
00242 q /= x;
00243 r__ /= x;
00244 L170:
00245 d__1 = sqrt(p * p + q * q + r__ * r__);
00246 s = d_sign(&d__1, &p);
00247 if (k == m) {
00248 goto L180;
00249 }
00250 h__[k + (k - 1) * h_dim1] = -s * x;
00251 goto L190;
00252 L180:
00253 if (l != m) {
00254 h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
00255 }
00256 L190:
00257 p += s;
00258 x = p / s;
00259 y = q / s;
00260 zz = r__ / s;
00261 q /= p;
00262 r__ /= p;
00263 if (notlas) {
00264 goto L225;
00265 }
00266
00267 i__2 = *n;
00268 for (j = k; j <= i__2; ++j) {
00269 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
00270 h__[k + j * h_dim1] -= p * x;
00271 h__[k + 1 + j * h_dim1] -= p * y;
00272
00273 }
00274
00275
00276 i__2 = en, i__3 = k + 3;
00277 j = min(i__2,i__3);
00278
00279 i__2 = j;
00280 for (i__ = 1; i__ <= i__2; ++i__) {
00281 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
00282 h__[i__ + k * h_dim1] -= p;
00283 h__[i__ + (k + 1) * h_dim1] -= p * q;
00284
00285 }
00286 goto L255;
00287 L225:
00288
00289 i__2 = *n;
00290 for (j = k; j <= i__2; ++j) {
00291 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
00292 k + 2 + j * h_dim1];
00293 h__[k + j * h_dim1] -= p * x;
00294 h__[k + 1 + j * h_dim1] -= p * y;
00295 h__[k + 2 + j * h_dim1] -= p * zz;
00296
00297 }
00298
00299
00300 i__2 = en, i__3 = k + 3;
00301 j = min(i__2,i__3);
00302
00303 i__2 = j;
00304 for (i__ = 1; i__ <= i__2; ++i__) {
00305 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] +
00306 zz * h__[i__ + (k + 2) * h_dim1];
00307 h__[i__ + k * h_dim1] -= p;
00308 h__[i__ + (k + 1) * h_dim1] -= p * q;
00309 h__[i__ + (k + 2) * h_dim1] -= p * r__;
00310
00311 }
00312 L255:
00313
00314 L260:
00315 ;
00316 }
00317
00318 goto L70;
00319
00320 L270:
00321 wr[en] = x + t;
00322 wi[en] = 0.;
00323 en = na;
00324 goto L60;
00325
00326 L280:
00327 p = (y - x) / 2.;
00328 q = p * p + w;
00329 zz = sqrt((abs(q)));
00330 x += t;
00331 if (q < 0.) {
00332 goto L320;
00333 }
00334
00335 zz = p + d_sign(&zz, &p);
00336 wr[na] = x + zz;
00337 wr[en] = wr[na];
00338 if (zz != 0.) {
00339 wr[en] = x - w / zz;
00340 }
00341 wi[na] = 0.;
00342 wi[en] = 0.;
00343 goto L330;
00344
00345 L320:
00346 wr[na] = x + p;
00347 wr[en] = x + p;
00348 wi[na] = zz;
00349 wi[en] = -zz;
00350 L330:
00351 en = enm2;
00352 goto L60;
00353
00354
00355 L1000:
00356 *ierr = en;
00357 L1001:
00358 return 0;
00359 }
00360