Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_bandr.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int bandr_(integer *nm, integer *n, integer *mb, doublereal *
00009 a, doublereal *d__, doublereal *e, doublereal *e2, logical *matz,
00010 doublereal *z__)
00011 {
00012
00013 integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5,
00014 i__6;
00015 doublereal d__1;
00016
00017
00018 double sqrt(doublereal);
00019
00020
00021 static doublereal dmin__;
00022 static integer maxl, maxr;
00023 static doublereal g;
00024 static integer j, k, l, r__;
00025 static doublereal u, b1, b2, c2, f1, f2;
00026 static integer i1, i2, j1, j2, m1, n2, r1;
00027 static doublereal s2;
00028 static integer kr, mr;
00029 static doublereal dminrt;
00030 static integer ugl;
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 z_dim1 = *nm;
00095 z_offset = z_dim1 + 1;
00096 z__ -= z_offset;
00097 --e2;
00098 --e;
00099 --d__;
00100 a_dim1 = *nm;
00101 a_offset = a_dim1 + 1;
00102 a -= a_offset;
00103
00104
00105 dmin__ = 5.4210108624275222e-20;
00106 dminrt = 2.3283064365386963e-10;
00107
00108 i__1 = *n;
00109 for (j = 1; j <= i__1; ++j) {
00110
00111 d__[j] = 1.;
00112 }
00113
00114 if (! (*matz)) {
00115 goto L60;
00116 }
00117
00118 i__1 = *n;
00119 for (j = 1; j <= i__1; ++j) {
00120
00121 i__2 = *n;
00122 for (k = 1; k <= i__2; ++k) {
00123
00124 z__[j + k * z_dim1] = 0.;
00125 }
00126
00127 z__[j + j * z_dim1] = 1.;
00128
00129 }
00130
00131 L60:
00132 m1 = *mb - 1;
00133 if ((i__1 = m1 - 1) < 0) {
00134 goto L900;
00135 } else if (i__1 == 0) {
00136 goto L800;
00137 } else {
00138 goto L70;
00139 }
00140 L70:
00141 n2 = *n - 2;
00142
00143 i__1 = n2;
00144 for (k = 1; k <= i__1; ++k) {
00145
00146 i__2 = m1, i__3 = *n - k;
00147 maxr = min(i__2,i__3);
00148
00149 i__2 = maxr;
00150 for (r1 = 2; r1 <= i__2; ++r1) {
00151 r__ = maxr + 2 - r1;
00152 kr = k + r__;
00153 mr = *mb - r__;
00154 g = a[kr + mr * a_dim1];
00155 a[kr - 1 + a_dim1] = a[kr - 1 + (mr + 1) * a_dim1];
00156 ugl = k;
00157
00158 i__3 = *n;
00159 i__4 = m1;
00160 for (j = kr; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
00161 j1 = j - 1;
00162 j2 = j1 - 1;
00163 if (g == 0.) {
00164 goto L600;
00165 }
00166 b1 = a[j1 + a_dim1] / g;
00167 b2 = b1 * d__[j1] / d__[j];
00168 s2 = 1. / (b1 * b2 + 1.);
00169 if (s2 >= .5) {
00170 goto L450;
00171 }
00172 b1 = g / a[j1 + a_dim1];
00173 b2 = b1 * d__[j] / d__[j1];
00174 c2 = 1. - s2;
00175 d__[j1] = c2 * d__[j1];
00176 d__[j] = c2 * d__[j];
00177 f1 = a[j + m1 * a_dim1] * 2.;
00178 f2 = b1 * a[j1 + *mb * a_dim1];
00179 a[j + m1 * a_dim1] = -b2 * (b1 * a[j + m1 * a_dim1] - a[j + *
00180 mb * a_dim1]) - f2 + a[j + m1 * a_dim1];
00181 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j + *mb * a_dim1] + f1) +
00182 a[j1 + *mb * a_dim1];
00183 a[j + *mb * a_dim1] = b1 * (f2 - f1) + a[j + *mb * a_dim1];
00184
00185 i__5 = j2;
00186 for (l = ugl; l <= i__5; ++l) {
00187 i2 = *mb - j + l;
00188 u = a[j1 + (i2 + 1) * a_dim1] + b2 * a[j + i2 * a_dim1];
00189 a[j + i2 * a_dim1] = -b1 * a[j1 + (i2 + 1) * a_dim1] + a[
00190 j + i2 * a_dim1];
00191 a[j1 + (i2 + 1) * a_dim1] = u;
00192
00193 }
00194
00195 ugl = j;
00196 a[j1 + a_dim1] += b2 * g;
00197 if (j == *n) {
00198 goto L350;
00199 }
00200
00201 i__5 = m1, i__6 = *n - j1;
00202 maxl = min(i__5,i__6);
00203
00204 i__5 = maxl;
00205 for (l = 2; l <= i__5; ++l) {
00206 i1 = j1 + l;
00207 i2 = *mb - l;
00208 u = a[i1 + i2 * a_dim1] + b2 * a[i1 + (i2 + 1) * a_dim1];
00209 a[i1 + (i2 + 1) * a_dim1] = -b1 * a[i1 + i2 * a_dim1] + a[
00210 i1 + (i2 + 1) * a_dim1];
00211 a[i1 + i2 * a_dim1] = u;
00212
00213 }
00214
00215 i1 = j + m1;
00216 if (i1 > *n) {
00217 goto L350;
00218 }
00219 g = b2 * a[i1 + a_dim1];
00220 L350:
00221 if (! (*matz)) {
00222 goto L500;
00223 }
00224
00225 i__5 = *n;
00226 for (l = 1; l <= i__5; ++l) {
00227 u = z__[l + j1 * z_dim1] + b2 * z__[l + j * z_dim1];
00228 z__[l + j * z_dim1] = -b1 * z__[l + j1 * z_dim1] + z__[l
00229 + j * z_dim1];
00230 z__[l + j1 * z_dim1] = u;
00231
00232 }
00233
00234 goto L500;
00235
00236 L450:
00237 u = d__[j1];
00238 d__[j1] = s2 * d__[j];
00239 d__[j] = s2 * u;
00240 f1 = a[j + m1 * a_dim1] * 2.;
00241 f2 = b1 * a[j + *mb * a_dim1];
00242 u = b1 * (f2 - f1) + a[j1 + *mb * a_dim1];
00243 a[j + m1 * a_dim1] = b2 * (b1 * a[j + m1 * a_dim1] - a[j1 + *
00244 mb * a_dim1]) + f2 - a[j + m1 * a_dim1];
00245 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j1 + *mb * a_dim1] + f1)
00246 + a[j + *mb * a_dim1];
00247 a[j + *mb * a_dim1] = u;
00248
00249 i__5 = j2;
00250 for (l = ugl; l <= i__5; ++l) {
00251 i2 = *mb - j + l;
00252 u = b2 * a[j1 + (i2 + 1) * a_dim1] + a[j + i2 * a_dim1];
00253 a[j + i2 * a_dim1] = -a[j1 + (i2 + 1) * a_dim1] + b1 * a[
00254 j + i2 * a_dim1];
00255 a[j1 + (i2 + 1) * a_dim1] = u;
00256
00257 }
00258
00259 ugl = j;
00260 a[j1 + a_dim1] = b2 * a[j1 + a_dim1] + g;
00261 if (j == *n) {
00262 goto L480;
00263 }
00264
00265 i__5 = m1, i__6 = *n - j1;
00266 maxl = min(i__5,i__6);
00267
00268 i__5 = maxl;
00269 for (l = 2; l <= i__5; ++l) {
00270 i1 = j1 + l;
00271 i2 = *mb - l;
00272 u = b2 * a[i1 + i2 * a_dim1] + a[i1 + (i2 + 1) * a_dim1];
00273 a[i1 + (i2 + 1) * a_dim1] = -a[i1 + i2 * a_dim1] + b1 * a[
00274 i1 + (i2 + 1) * a_dim1];
00275 a[i1 + i2 * a_dim1] = u;
00276
00277 }
00278
00279 i1 = j + m1;
00280 if (i1 > *n) {
00281 goto L480;
00282 }
00283 g = a[i1 + a_dim1];
00284 a[i1 + a_dim1] = b1 * a[i1 + a_dim1];
00285 L480:
00286 if (! (*matz)) {
00287 goto L500;
00288 }
00289
00290 i__5 = *n;
00291 for (l = 1; l <= i__5; ++l) {
00292 u = b2 * z__[l + j1 * z_dim1] + z__[l + j * z_dim1];
00293 z__[l + j * z_dim1] = -z__[l + j1 * z_dim1] + b1 * z__[l
00294 + j * z_dim1];
00295 z__[l + j1 * z_dim1] = u;
00296
00297 }
00298
00299 L500:
00300 ;
00301 }
00302
00303 L600:
00304 ;
00305 }
00306
00307 if (k % 64 != 0) {
00308 goto L700;
00309 }
00310
00311 i__2 = *n;
00312 for (j = k; j <= i__2; ++j) {
00313 if (d__[j] >= dmin__) {
00314 goto L650;
00315 }
00316
00317 i__4 = 1, i__3 = *mb + 1 - j;
00318 maxl = max(i__4,i__3);
00319
00320 i__4 = m1;
00321 for (l = maxl; l <= i__4; ++l) {
00322
00323 a[j + l * a_dim1] = dminrt * a[j + l * a_dim1];
00324 }
00325
00326 if (j == *n) {
00327 goto L630;
00328 }
00329
00330 i__4 = m1, i__3 = *n - j;
00331 maxl = min(i__4,i__3);
00332
00333 i__4 = maxl;
00334 for (l = 1; l <= i__4; ++l) {
00335 i1 = j + l;
00336 i2 = *mb - l;
00337 a[i1 + i2 * a_dim1] = dminrt * a[i1 + i2 * a_dim1];
00338
00339 }
00340
00341 L630:
00342 if (! (*matz)) {
00343 goto L645;
00344 }
00345
00346 i__4 = *n;
00347 for (l = 1; l <= i__4; ++l) {
00348
00349 z__[l + j * z_dim1] = dminrt * z__[l + j * z_dim1];
00350 }
00351
00352 L645:
00353 a[j + *mb * a_dim1] = dmin__ * a[j + *mb * a_dim1];
00354 d__[j] /= dmin__;
00355 L650:
00356 ;
00357 }
00358
00359 L700:
00360 ;
00361 }
00362
00363 L800:
00364 i__1 = *n;
00365 for (j = 2; j <= i__1; ++j) {
00366
00367 e[j] = sqrt(d__[j]);
00368 }
00369
00370 if (! (*matz)) {
00371 goto L840;
00372 }
00373
00374 i__1 = *n;
00375 for (j = 1; j <= i__1; ++j) {
00376
00377 i__2 = *n;
00378 for (k = 2; k <= i__2; ++k) {
00379
00380 z__[j + k * z_dim1] = e[k] * z__[j + k * z_dim1];
00381 }
00382
00383
00384 }
00385
00386 L840:
00387 u = 1.;
00388
00389 i__1 = *n;
00390 for (j = 2; j <= i__1; ++j) {
00391 a[j + m1 * a_dim1] = u * e[j] * a[j + m1 * a_dim1];
00392 u = e[j];
00393
00394 d__1 = a[j + m1 * a_dim1];
00395 e2[j] = d__1 * d__1;
00396 a[j + *mb * a_dim1] = d__[j] * a[j + *mb * a_dim1];
00397 d__[j] = a[j + *mb * a_dim1];
00398 e[j] = a[j + m1 * a_dim1];
00399
00400 }
00401
00402 d__[1] = a[*mb * a_dim1 + 1];
00403 e[1] = 0.;
00404 e2[1] = 0.;
00405 goto L1001;
00406
00407 L900:
00408 i__1 = *n;
00409 for (j = 1; j <= i__1; ++j) {
00410 d__[j] = a[j + *mb * a_dim1];
00411 e[j] = 0.;
00412 e2[j] = 0.;
00413
00414 }
00415
00416 L1001:
00417 return 0;
00418 }
00419