Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_bandv.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int bandv_(integer *nm, integer *n, integer *mbw, doublereal
00009 *a, doublereal *e21, integer *m, doublereal *w, doublereal *z__,
00010 integer *ierr, integer *nv, doublereal *rv, doublereal *rv6)
00011 {
00012
00013 integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
00014 doublereal d__1;
00015
00016
00017 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00018
00019
00020 static integer maxj, maxk;
00021 static doublereal norm;
00022 static integer i__, j, k, r__;
00023 static doublereal u, v, order;
00024 static integer group, m1;
00025 static doublereal x0, x1;
00026 static integer mb, m21, ii, ij, jj, kj;
00027 static doublereal uk, xu;
00028 extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal
00029 *);
00030 static integer ij1, kj1, its;
00031 static doublereal eps2, eps3, eps4;
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
00127
00128
00129
00130
00131
00132
00133
00134
00135 --rv6;
00136 a_dim1 = *nm;
00137 a_offset = a_dim1 + 1;
00138 a -= a_offset;
00139 z_dim1 = *nm;
00140 z_offset = z_dim1 + 1;
00141 z__ -= z_offset;
00142 --w;
00143 --rv;
00144
00145
00146 *ierr = 0;
00147 if (*m == 0) {
00148 goto L1001;
00149 }
00150 mb = *mbw;
00151 if (*e21 < 0.) {
00152 mb = (*mbw + 1) / 2;
00153 }
00154 m1 = mb - 1;
00155 m21 = m1 + mb;
00156 order = 1. - abs(*e21);
00157
00158 i__1 = *m;
00159 for (r__ = 1; r__ <= i__1; ++r__) {
00160 its = 1;
00161 x1 = w[r__];
00162 if (r__ != 1) {
00163 goto L100;
00164 }
00165
00166 norm = 0.;
00167
00168 i__2 = mb;
00169 for (j = 1; j <= i__2; ++j) {
00170 jj = mb + 1 - j;
00171 kj = jj + m1;
00172 ij = 1;
00173 v = 0.;
00174
00175 i__3 = *n;
00176 for (i__ = jj; i__ <= i__3; ++i__) {
00177 v += (d__1 = a[i__ + j * a_dim1], abs(d__1));
00178 if (*e21 >= 0.) {
00179 goto L40;
00180 }
00181 v += (d__1 = a[ij + kj * a_dim1], abs(d__1));
00182 ++ij;
00183 L40:
00184 ;
00185 }
00186
00187 norm = max(norm,v);
00188
00189 }
00190
00191 if (*e21 < 0.) {
00192 norm *= .5;
00193 }
00194
00195
00196
00197
00198
00199 if (norm == 0.) {
00200 norm = 1.;
00201 }
00202 eps2 = norm * .001 * abs(order);
00203 eps3 = epslon_(&norm);
00204 uk = (doublereal) (*n);
00205 uk = sqrt(uk);
00206 eps4 = uk * eps3;
00207 L80:
00208 group = 0;
00209 goto L120;
00210
00211 L100:
00212 if ((d__1 = x1 - x0, abs(d__1)) >= eps2) {
00213 goto L80;
00214 }
00215 ++group;
00216 if (order * (x1 - x0) <= 0.) {
00217 x1 = x0 + order * eps3;
00218 }
00219
00220
00221 L120:
00222 i__2 = *n;
00223 for (i__ = 1; i__ <= i__2; ++i__) {
00224
00225 i__3 = 0, i__4 = i__ - m1;
00226 ij = i__ + min(i__3,i__4) * *n;
00227 kj = ij + mb * *n;
00228 ij1 = kj + m1 * *n;
00229 if (m1 == 0) {
00230 goto L180;
00231 }
00232
00233 i__3 = m1;
00234 for (j = 1; j <= i__3; ++j) {
00235 if (ij > m1) {
00236 goto L125;
00237 }
00238 if (ij > 0) {
00239 goto L130;
00240 }
00241 rv[ij1] = 0.;
00242 ij1 += *n;
00243 goto L130;
00244 L125:
00245 rv[ij] = a[i__ + j * a_dim1];
00246 L130:
00247 ij += *n;
00248 ii = i__ + j;
00249 if (ii > *n) {
00250 goto L150;
00251 }
00252 jj = mb - j;
00253 if (*e21 >= 0.) {
00254 goto L140;
00255 }
00256 ii = i__;
00257 jj = mb + j;
00258 L140:
00259 rv[kj] = a[ii + jj * a_dim1];
00260 kj += *n;
00261 L150:
00262 ;
00263 }
00264
00265 L180:
00266 rv[ij] = a[i__ + mb * a_dim1] - x1;
00267 rv6[i__] = eps4;
00268 if (order == 0.) {
00269 rv6[i__] = z__[i__ + r__ * z_dim1];
00270 }
00271
00272 }
00273
00274 if (m1 == 0) {
00275 goto L600;
00276 }
00277
00278 i__2 = *n;
00279 for (i__ = 1; i__ <= i__2; ++i__) {
00280 ii = i__ + 1;
00281
00282 i__3 = i__ + m1 - 1;
00283 maxk = min(i__3,*n);
00284
00285 i__3 = *n - i__, i__4 = m21 - 2;
00286 maxj = min(i__3,i__4) * *n;
00287
00288 i__3 = maxk;
00289 for (k = i__; k <= i__3; ++k) {
00290 kj1 = k;
00291 j = kj1 + *n;
00292 jj = j + maxj;
00293
00294 i__4 = jj;
00295 i__5 = *n;
00296 for (kj = j; i__5 < 0 ? kj >= i__4 : kj <= i__4; kj += i__5) {
00297 rv[kj1] = rv[kj];
00298 kj1 = kj;
00299
00300 }
00301
00302 rv[kj1] = 0.;
00303
00304 }
00305
00306 if (i__ == *n) {
00307 goto L580;
00308 }
00309 u = 0.;
00310
00311 i__3 = i__ + m1;
00312 maxk = min(i__3,*n);
00313
00314 i__3 = *n - ii, i__5 = m21 - 2;
00315 maxj = min(i__3,i__5) * *n;
00316
00317 i__3 = maxk;
00318 for (j = i__; j <= i__3; ++j) {
00319 if ((d__1 = rv[j], abs(d__1)) < abs(u)) {
00320 goto L450;
00321 }
00322 u = rv[j];
00323 k = j;
00324 L450:
00325 ;
00326 }
00327
00328 j = i__ + *n;
00329 jj = j + maxj;
00330 if (k == i__) {
00331 goto L520;
00332 }
00333 kj = k;
00334
00335 i__3 = jj;
00336 i__5 = *n;
00337 for (ij = i__; i__5 < 0 ? ij >= i__3 : ij <= i__3; ij += i__5) {
00338 v = rv[ij];
00339 rv[ij] = rv[kj];
00340 rv[kj] = v;
00341 kj += *n;
00342
00343 }
00344
00345 if (order != 0.) {
00346 goto L520;
00347 }
00348 v = rv6[i__];
00349 rv6[i__] = rv6[k];
00350 rv6[k] = v;
00351 L520:
00352 if (u == 0.) {
00353 goto L580;
00354 }
00355
00356 i__5 = maxk;
00357 for (k = ii; k <= i__5; ++k) {
00358 v = rv[k] / u;
00359 kj = k;
00360
00361 i__3 = jj;
00362 i__4 = *n;
00363 for (ij = j; i__4 < 0 ? ij >= i__3 : ij <= i__3; ij += i__4) {
00364 kj += *n;
00365 rv[kj] -= v * rv[ij];
00366
00367 }
00368
00369 if (order == 0.) {
00370 rv6[k] -= v * rv6[i__];
00371 }
00372
00373 }
00374
00375 L580:
00376 ;
00377 }
00378
00379
00380 L600:
00381 i__2 = *n;
00382 for (ii = 1; ii <= i__2; ++ii) {
00383 i__ = *n + 1 - ii;
00384 maxj = min(ii,m21);
00385 if (maxj == 1) {
00386 goto L620;
00387 }
00388 ij1 = i__;
00389 j = ij1 + *n;
00390 jj = j + (maxj - 2) * *n;
00391
00392 i__5 = jj;
00393 i__4 = *n;
00394 for (ij = j; i__4 < 0 ? ij >= i__5 : ij <= i__5; ij += i__4) {
00395 ++ij1;
00396 rv6[i__] -= rv[ij] * rv6[ij1];
00397
00398 }
00399
00400 L620:
00401 v = rv[i__];
00402 if (abs(v) >= eps3) {
00403 goto L625;
00404 }
00405
00406
00407 if (order == 0.) {
00408 *ierr = -r__;
00409 }
00410 v = d_sign(&eps3, &v);
00411 L625:
00412 rv6[i__] /= v;
00413
00414 }
00415
00416 xu = 1.;
00417 if (order == 0.) {
00418 goto L870;
00419 }
00420
00421
00422 if (group == 0) {
00423 goto L700;
00424 }
00425
00426 i__2 = group;
00427 for (jj = 1; jj <= i__2; ++jj) {
00428 j = r__ - group - 1 + jj;
00429 xu = 0.;
00430
00431 i__4 = *n;
00432 for (i__ = 1; i__ <= i__4; ++i__) {
00433
00434 xu += rv6[i__] * z__[i__ + j * z_dim1];
00435 }
00436
00437 i__4 = *n;
00438 for (i__ = 1; i__ <= i__4; ++i__) {
00439
00440 rv6[i__] -= xu * z__[i__ + j * z_dim1];
00441 }
00442
00443
00444 }
00445
00446 L700:
00447 norm = 0.;
00448
00449 i__2 = *n;
00450 for (i__ = 1; i__ <= i__2; ++i__) {
00451
00452 norm += (d__1 = rv6[i__], abs(d__1));
00453 }
00454
00455 if (norm >= .1) {
00456 goto L840;
00457 }
00458
00459
00460 if (its >= *n) {
00461 goto L830;
00462 }
00463 ++its;
00464 xu = eps4 / (uk + 1.);
00465 rv6[1] = eps4;
00466
00467 i__2 = *n;
00468 for (i__ = 2; i__ <= i__2; ++i__) {
00469
00470 rv6[i__] = xu;
00471 }
00472
00473 rv6[its] -= eps4 * uk;
00474 goto L600;
00475
00476 L830:
00477 *ierr = -r__;
00478 xu = 0.;
00479 goto L870;
00480
00481
00482 L840:
00483 u = 0.;
00484
00485 i__2 = *n;
00486 for (i__ = 1; i__ <= i__2; ++i__) {
00487
00488 u = pythag_(&u, &rv6[i__]);
00489 }
00490
00491 xu = 1. / u;
00492
00493 L870:
00494 i__2 = *n;
00495 for (i__ = 1; i__ <= i__2; ++i__) {
00496
00497 z__[i__ + r__ * z_dim1] = rv6[i__] * xu;
00498 }
00499
00500 x0 = x1;
00501
00502 }
00503
00504 L1001:
00505 return 0;
00506 }
00507