Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_tsturm.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008
00009
00010 static doublereal c_b26 = 1.;
00011
00012 int tsturm_(integer *nm, integer *n, doublereal *eps1,
00013 doublereal *d__, doublereal *e, doublereal *e2, doublereal *lb,
00014 doublereal *ub, integer *mm, integer *m, doublereal *w, doublereal *
00015 z__, integer *ierr, doublereal *rv1, doublereal *rv2, doublereal *rv3,
00016 doublereal *rv4, doublereal *rv5, doublereal *rv6)
00017 {
00018
00019 integer z_dim1, z_offset, i__1, i__2, i__3;
00020 doublereal d__1, d__2, d__3, d__4;
00021
00022
00023 double sqrt(doublereal);
00024
00025
00026 static doublereal norm;
00027 static integer i__, j, k, p, q, r__, s;
00028 static doublereal u, v;
00029 static integer group, m1, m2;
00030 static doublereal t1, t2, x0, x1;
00031 static integer ii, jj, ip;
00032 static doublereal uk, xu;
00033 extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal
00034 *);
00035 static integer isturm, its;
00036 static doublereal eps2, eps3, eps4, tst1, tst2;
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 --rv6;
00129 --rv5;
00130 --rv4;
00131 --rv3;
00132 --rv2;
00133 --rv1;
00134 --e2;
00135 --e;
00136 --d__;
00137 z_dim1 = *nm;
00138 z_offset = z_dim1 + 1;
00139 z__ -= z_offset;
00140 --w;
00141
00142
00143 *ierr = 0;
00144 t1 = *lb;
00145 t2 = *ub;
00146
00147 i__1 = *n;
00148 for (i__ = 1; i__ <= i__1; ++i__) {
00149 if (i__ == 1) {
00150 goto L20;
00151 }
00152 tst1 = (d__1 = d__[i__], abs(d__1)) + (d__2 = d__[i__ - 1], abs(d__2))
00153 ;
00154 tst2 = tst1 + (d__1 = e[i__], abs(d__1));
00155 if (tst2 > tst1) {
00156 goto L40;
00157 }
00158 L20:
00159 e2[i__] = 0.;
00160 L40:
00161 ;
00162 }
00163
00164
00165 p = 1;
00166 q = *n;
00167 x1 = *ub;
00168 isturm = 1;
00169 goto L320;
00170 L60:
00171 *m = s;
00172 x1 = *lb;
00173 isturm = 2;
00174 goto L320;
00175 L80:
00176 *m -= s;
00177 if (*m > *mm) {
00178 goto L980;
00179 }
00180 q = 0;
00181 r__ = 0;
00182
00183
00184 L100:
00185 if (r__ == *m) {
00186 goto L1001;
00187 }
00188 p = q + 1;
00189 xu = d__[p];
00190 x0 = d__[p];
00191 u = 0.;
00192
00193 i__1 = *n;
00194 for (q = p; q <= i__1; ++q) {
00195 x1 = u;
00196 u = 0.;
00197 v = 0.;
00198 if (q == *n) {
00199 goto L110;
00200 }
00201 u = (d__1 = e[q + 1], abs(d__1));
00202 v = e2[q + 1];
00203 L110:
00204
00205 d__1 = d__[q] - (x1 + u);
00206 xu = min(d__1,xu);
00207
00208 d__1 = d__[q] + (x1 + u);
00209 x0 = max(d__1,x0);
00210 if (v == 0.) {
00211 goto L140;
00212 }
00213
00214 }
00215
00216 L140:
00217
00218 d__2 = abs(xu), d__3 = abs(x0);
00219 d__1 = max(d__2,d__3);
00220 x1 = epslon_(&d__1);
00221 if (*eps1 <= 0.) {
00222 *eps1 = -x1;
00223 }
00224 if (p != q) {
00225 goto L180;
00226 }
00227
00228 if (t1 > d__[p] || d__[p] >= t2) {
00229 goto L940;
00230 }
00231 ++r__;
00232
00233 i__1 = *n;
00234 for (i__ = 1; i__ <= i__1; ++i__) {
00235
00236 z__[i__ + r__ * z_dim1] = 0.;
00237 }
00238
00239 w[r__] = d__[p];
00240 z__[p + r__ * z_dim1] = 1.;
00241 goto L940;
00242 L180:
00243 u = (doublereal) (q - p + 1);
00244 x1 = u * x1;
00245
00246 d__1 = t1, d__2 = xu - x1;
00247 *lb = max(d__1,d__2);
00248
00249 d__1 = t2, d__2 = x0 + x1;
00250 *ub = min(d__1,d__2);
00251 x1 = *lb;
00252 isturm = 3;
00253 goto L320;
00254 L200:
00255 m1 = s + 1;
00256 x1 = *ub;
00257 isturm = 4;
00258 goto L320;
00259 L220:
00260 m2 = s;
00261 if (m1 > m2) {
00262 goto L940;
00263 }
00264
00265 x0 = *ub;
00266 isturm = 5;
00267
00268 i__1 = m2;
00269 for (i__ = m1; i__ <= i__1; ++i__) {
00270 rv5[i__] = *ub;
00271 rv4[i__] = *lb;
00272
00273 }
00274
00275
00276
00277
00278 k = m2;
00279 L250:
00280 xu = *lb;
00281
00282 i__1 = k;
00283 for (ii = m1; ii <= i__1; ++ii) {
00284 i__ = m1 + k - ii;
00285 if (xu >= rv4[i__]) {
00286 goto L260;
00287 }
00288 xu = rv4[i__];
00289 goto L280;
00290 L260:
00291 ;
00292 }
00293
00294 L280:
00295 if (x0 > rv5[k]) {
00296 x0 = rv5[k];
00297 }
00298
00299 L300:
00300 x1 = (xu + x0) * .5;
00301 if (x0 - xu <= abs(*eps1)) {
00302 goto L420;
00303 }
00304 tst1 = (abs(xu) + abs(x0)) * 2.;
00305 tst2 = tst1 + (x0 - xu);
00306 if (tst2 == tst1) {
00307 goto L420;
00308 }
00309
00310 L320:
00311 s = p - 1;
00312 u = 1.;
00313
00314 i__1 = q;
00315 for (i__ = p; i__ <= i__1; ++i__) {
00316 if (u != 0.) {
00317 goto L325;
00318 }
00319 v = (d__1 = e[i__], abs(d__1)) / epslon_(&c_b26);
00320 if (e2[i__] == 0.) {
00321 v = 0.;
00322 }
00323 goto L330;
00324 L325:
00325 v = e2[i__] / u;
00326 L330:
00327 u = d__[i__] - x1 - v;
00328 if (u < 0.) {
00329 ++s;
00330 }
00331
00332 }
00333
00334 switch (isturm) {
00335 case 1: goto L60;
00336 case 2: goto L80;
00337 case 3: goto L200;
00338 case 4: goto L220;
00339 case 5: goto L360;
00340 }
00341
00342 L360:
00343 if (s >= k) {
00344 goto L400;
00345 }
00346 xu = x1;
00347 if (s >= m1) {
00348 goto L380;
00349 }
00350 rv4[m1] = x1;
00351 goto L300;
00352 L380:
00353 rv4[s + 1] = x1;
00354 if (rv5[s] > x1) {
00355 rv5[s] = x1;
00356 }
00357 goto L300;
00358 L400:
00359 x0 = x1;
00360 goto L300;
00361
00362 L420:
00363 rv5[k] = x1;
00364 --k;
00365 if (k >= m1) {
00366 goto L250;
00367 }
00368
00369 norm = (d__1 = d__[p], abs(d__1));
00370 ip = p + 1;
00371
00372 i__1 = q;
00373 for (i__ = ip; i__ <= i__1; ++i__) {
00374
00375
00376 d__3 = norm, d__4 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[i__],
00377 abs(d__2));
00378 norm = max(d__3,d__4);
00379 }
00380
00381
00382
00383
00384 eps2 = norm * .001;
00385 eps3 = epslon_(&norm);
00386 uk = (doublereal) (q - p + 1);
00387 eps4 = uk * eps3;
00388 uk = eps4 / sqrt(uk);
00389 group = 0;
00390 s = p;
00391
00392 i__1 = m2;
00393 for (k = m1; k <= i__1; ++k) {
00394 ++r__;
00395 its = 1;
00396 w[r__] = rv5[k];
00397 x1 = rv5[k];
00398
00399 if (k == m1) {
00400 goto L520;
00401 }
00402 if (x1 - x0 >= eps2) {
00403 group = -1;
00404 }
00405 ++group;
00406 if (x1 <= x0) {
00407 x1 = x0 + eps3;
00408 }
00409
00410
00411 L520:
00412 v = 0.;
00413
00414 i__2 = q;
00415 for (i__ = p; i__ <= i__2; ++i__) {
00416 rv6[i__] = uk;
00417 if (i__ == p) {
00418 goto L560;
00419 }
00420 if ((d__1 = e[i__], abs(d__1)) < abs(u)) {
00421 goto L540;
00422 }
00423 xu = u / e[i__];
00424 rv4[i__] = xu;
00425 rv1[i__ - 1] = e[i__];
00426 rv2[i__ - 1] = d__[i__] - x1;
00427 rv3[i__ - 1] = 0.;
00428 if (i__ != q) {
00429 rv3[i__ - 1] = e[i__ + 1];
00430 }
00431 u = v - xu * rv2[i__ - 1];
00432 v = -xu * rv3[i__ - 1];
00433 goto L580;
00434 L540:
00435 xu = e[i__] / u;
00436 rv4[i__] = xu;
00437 rv1[i__ - 1] = u;
00438 rv2[i__ - 1] = v;
00439 rv3[i__ - 1] = 0.;
00440 L560:
00441 u = d__[i__] - x1 - xu * v;
00442 if (i__ != q) {
00443 v = e[i__ + 1];
00444 }
00445 L580:
00446 ;
00447 }
00448
00449 if (u == 0.) {
00450 u = eps3;
00451 }
00452 rv1[q] = u;
00453 rv2[q] = 0.;
00454 rv3[q] = 0.;
00455
00456
00457 L600:
00458 i__2 = q;
00459 for (ii = p; ii <= i__2; ++ii) {
00460 i__ = p + q - ii;
00461 rv6[i__] = (rv6[i__] - u * rv2[i__] - v * rv3[i__]) / rv1[i__];
00462 v = u;
00463 u = rv6[i__];
00464
00465 }
00466
00467
00468 if (group == 0) {
00469 goto L700;
00470 }
00471
00472 i__2 = group;
00473 for (jj = 1; jj <= i__2; ++jj) {
00474 j = r__ - group - 1 + jj;
00475 xu = 0.;
00476
00477 i__3 = q;
00478 for (i__ = p; i__ <= i__3; ++i__) {
00479
00480 xu += rv6[i__] * z__[i__ + j * z_dim1];
00481 }
00482
00483 i__3 = q;
00484 for (i__ = p; i__ <= i__3; ++i__) {
00485
00486 rv6[i__] -= xu * z__[i__ + j * z_dim1];
00487 }
00488
00489
00490 }
00491
00492 L700:
00493 norm = 0.;
00494
00495 i__2 = q;
00496 for (i__ = p; i__ <= i__2; ++i__) {
00497
00498 norm += (d__1 = rv6[i__], abs(d__1));
00499 }
00500
00501 if (norm >= 1.) {
00502 goto L840;
00503 }
00504
00505 if (its == 5) {
00506 goto L960;
00507 }
00508 if (norm != 0.) {
00509 goto L740;
00510 }
00511 rv6[s] = eps4;
00512 ++s;
00513 if (s > q) {
00514 s = p;
00515 }
00516 goto L780;
00517 L740:
00518 xu = eps4 / norm;
00519
00520 i__2 = q;
00521 for (i__ = p; i__ <= i__2; ++i__) {
00522
00523 rv6[i__] *= xu;
00524 }
00525
00526
00527 L780:
00528 i__2 = q;
00529 for (i__ = ip; i__ <= i__2; ++i__) {
00530 u = rv6[i__];
00531
00532
00533
00534 if (rv1[i__ - 1] != e[i__]) {
00535 goto L800;
00536 }
00537 u = rv6[i__ - 1];
00538 rv6[i__ - 1] = rv6[i__];
00539 L800:
00540 rv6[i__] = u - rv4[i__] * rv6[i__ - 1];
00541
00542 }
00543
00544 ++its;
00545 goto L600;
00546
00547
00548 L840:
00549 u = 0.;
00550
00551 i__2 = q;
00552 for (i__ = p; i__ <= i__2; ++i__) {
00553
00554 u = pythag_(&u, &rv6[i__]);
00555 }
00556
00557 xu = 1. / u;
00558
00559 i__2 = *n;
00560 for (i__ = 1; i__ <= i__2; ++i__) {
00561
00562 z__[i__ + r__ * z_dim1] = 0.;
00563 }
00564
00565 i__2 = q;
00566 for (i__ = p; i__ <= i__2; ++i__) {
00567
00568 z__[i__ + r__ * z_dim1] = rv6[i__] * xu;
00569 }
00570
00571 x0 = x1;
00572
00573 }
00574
00575 L940:
00576 if (q < *n) {
00577 goto L100;
00578 }
00579 goto L1001;
00580
00581 L960:
00582 *ierr = (*n << 2) + r__;
00583 goto L1001;
00584
00585
00586 L980:
00587 *ierr = *n * 3 + 1;
00588 L1001:
00589 *lb = t1;
00590 *ub = t2;
00591 return 0;
00592 }
00593