Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_bqr.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008
00009
00010 static doublereal c_b8 = 1.;
00011
00012 int bqr_(integer *nm, integer *n, integer *mb, doublereal *a,
00013 doublereal *t, doublereal *r__, integer *ierr, integer *nv,
00014 doublereal *rv)
00015 {
00016
00017 integer a_dim1, a_offset, i__1, i__2, i__3;
00018 doublereal d__1;
00019
00020
00021 double d_sign(doublereal *, doublereal *), sqrt(doublereal);
00022
00023
00024 static doublereal f, g;
00025 static integer i__, j, k, l, m;
00026 static doublereal q, s, scale;
00027 static integer imult, m1, m2, m3, m4, m21, m31, ii, ik, jk, kj, jm, kk,
00028 km, ll, mk, mn, ni, mz;
00029 extern doublereal pythag_(doublereal *, doublereal *);
00030 static integer kj1, its;
00031 static doublereal tst1, tst2;
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
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126 a_dim1 = *nm;
00127 a_offset = a_dim1 + 1;
00128 a -= a_offset;
00129 --rv;
00130
00131
00132 *ierr = 0;
00133 m1 = min(*mb,*n);
00134 m = m1 - 1;
00135 m2 = m + m;
00136 m21 = m2 + 1;
00137 m3 = m21 + m;
00138 m31 = m3 + 1;
00139 m4 = m31 + m2;
00140 mn = m + *n;
00141 mz = *mb - m1;
00142 its = 0;
00143
00144 L40:
00145 g = a[*n + *mb * a_dim1];
00146 if (m == 0) {
00147 goto L360;
00148 }
00149 f = 0.;
00150
00151 i__1 = m;
00152 for (k = 1; k <= i__1; ++k) {
00153 mk = k + mz;
00154 f += (d__1 = a[*n + mk * a_dim1], abs(d__1));
00155
00156 }
00157
00158 if (its == 0 && f > *r__) {
00159 *r__ = f;
00160 }
00161 tst1 = *r__;
00162 tst2 = tst1 + f;
00163 if (tst2 <= tst1) {
00164 goto L360;
00165 }
00166 if (its == 30) {
00167 goto L1000;
00168 }
00169 ++its;
00170
00171 if (f > *r__ * .25 && its < 5) {
00172 goto L90;
00173 }
00174 f = a[*n + (*mb - 1) * a_dim1];
00175 if (f == 0.) {
00176 goto L70;
00177 }
00178 q = (a[*n - 1 + *mb * a_dim1] - g) / (f * 2.);
00179 s = pythag_(&q, &c_b8);
00180 g -= f / (q + d_sign(&s, &q));
00181 L70:
00182 *t += g;
00183
00184 i__1 = *n;
00185 for (i__ = 1; i__ <= i__1; ++i__) {
00186
00187 a[i__ + *mb * a_dim1] -= g;
00188 }
00189
00190 L90:
00191 i__1 = m4;
00192 for (k = m31; k <= i__1; ++k) {
00193
00194 rv[k] = 0.;
00195 }
00196
00197 i__1 = mn;
00198 for (ii = 1; ii <= i__1; ++ii) {
00199 i__ = ii - m;
00200 ni = *n - ii;
00201 if (ni < 0) {
00202 goto L230;
00203 }
00204
00205
00206 i__2 = 1, i__3 = 2 - i__;
00207 l = max(i__2,i__3);
00208
00209 i__2 = m3;
00210 for (k = 1; k <= i__2; ++k) {
00211
00212 rv[k] = 0.;
00213 }
00214
00215 i__2 = m1;
00216 for (k = l; k <= i__2; ++k) {
00217 km = k + m;
00218 mk = k + mz;
00219 rv[km] = a[ii + mk * a_dim1];
00220
00221 }
00222
00223 ll = min(m,ni);
00224 if (ll == 0) {
00225 goto L135;
00226 }
00227
00228 i__2 = ll;
00229 for (k = 1; k <= i__2; ++k) {
00230 km = k + m21;
00231 ik = ii + k;
00232 mk = *mb - k;
00233 rv[km] = a[ik + mk * a_dim1];
00234
00235 }
00236
00237
00238 L135:
00239 ll = m2;
00240 imult = 0;
00241
00242 L140:
00243 kj = m4 - m1;
00244
00245 i__2 = ll;
00246 for (j = 1; j <= i__2; ++j) {
00247 kj += m1;
00248 jm = j + m3;
00249 if (rv[jm] == 0.) {
00250 goto L170;
00251 }
00252 f = 0.;
00253
00254 i__3 = m1;
00255 for (k = 1; k <= i__3; ++k) {
00256 ++kj;
00257 jk = j + k - 1;
00258 f += rv[kj] * rv[jk];
00259
00260 }
00261
00262 f /= rv[jm];
00263 kj -= m1;
00264
00265 i__3 = m1;
00266 for (k = 1; k <= i__3; ++k) {
00267 ++kj;
00268 jk = j + k - 1;
00269 rv[jk] -= rv[kj] * f;
00270
00271 }
00272
00273 kj -= m1;
00274 L170:
00275 ;
00276 }
00277
00278 if (imult != 0) {
00279 goto L280;
00280 }
00281
00282 f = rv[m21];
00283 s = 0.;
00284 rv[m4] = 0.;
00285 scale = 0.;
00286
00287 i__2 = m3;
00288 for (k = m21; k <= i__2; ++k) {
00289
00290 scale += (d__1 = rv[k], abs(d__1));
00291 }
00292
00293 if (scale == 0.) {
00294 goto L210;
00295 }
00296
00297 i__2 = m3;
00298 for (k = m21; k <= i__2; ++k) {
00299
00300
00301 d__1 = rv[k] / scale;
00302 s += d__1 * d__1;
00303 }
00304
00305 s = scale * scale * s;
00306 d__1 = sqrt(s);
00307 g = -d_sign(&d__1, &f);
00308 rv[m21] = g;
00309 rv[m4] = s - f * g;
00310 kj = m4 + m2 * m1 + 1;
00311 rv[kj] = f - g;
00312
00313 i__2 = m1;
00314 for (k = 2; k <= i__2; ++k) {
00315 ++kj;
00316 km = k + m2;
00317 rv[kj] = rv[km];
00318
00319 }
00320
00321 L210:
00322 i__2 = m1;
00323 for (k = l; k <= i__2; ++k) {
00324 km = k + m;
00325 mk = k + mz;
00326 a[ii + mk * a_dim1] = rv[km];
00327
00328 }
00329
00330 L230:
00331
00332 i__2 = 1, i__3 = m1 + 1 - i__;
00333 l = max(i__2,i__3);
00334 if (i__ <= 0) {
00335 goto L300;
00336 }
00337
00338 i__2 = m21;
00339 for (k = 1; k <= i__2; ++k) {
00340
00341 rv[k] = 0.;
00342 }
00343
00344
00345 i__2 = m1, i__3 = ni + m1;
00346 ll = min(i__2,i__3);
00347
00348 i__2 = ll;
00349 for (kk = 1; kk <= i__2; ++kk) {
00350 k = kk - 1;
00351 km = k + m1;
00352 ik = i__ + k;
00353 mk = *mb - k;
00354 rv[km] = a[ik + mk * a_dim1];
00355
00356 }
00357
00358
00359 ll = m1;
00360 imult = 1;
00361 goto L140;
00362
00363 L280:
00364 i__2 = m1;
00365 for (k = l; k <= i__2; ++k) {
00366 mk = k + mz;
00367 a[i__ + mk * a_dim1] = rv[k];
00368
00369 }
00370
00371 L300:
00372 if (l > 1) {
00373 --l;
00374 }
00375 kj1 = m4 + l * m1;
00376
00377 i__2 = m2;
00378 for (j = l; j <= i__2; ++j) {
00379 jm = j + m3;
00380 rv[jm] = rv[jm + 1];
00381
00382 i__3 = m1;
00383 for (k = 1; k <= i__3; ++k) {
00384 ++kj1;
00385 kj = kj1 - m1;
00386 rv[kj] = rv[kj1];
00387
00388 }
00389 }
00390
00391
00392 }
00393
00394 goto L40;
00395
00396 L360:
00397 *t += g;
00398
00399 i__1 = *n;
00400 for (i__ = 1; i__ <= i__1; ++i__) {
00401
00402 a[i__ + *mb * a_dim1] -= g;
00403 }
00404
00405 i__1 = m1;
00406 for (k = 1; k <= i__1; ++k) {
00407 mk = k + mz;
00408 a[*n + mk * a_dim1] = 0.;
00409
00410 }
00411
00412 goto L1001;
00413
00414
00415 L1000:
00416 *ierr = *n;
00417 L1001:
00418 return 0;
00419 }
00420