Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_bisect.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 bisect_(integer *n, doublereal *eps1, doublereal *d__,
00013 doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub,
00014 integer *mm, integer *m, doublereal *w, integer *ind, integer *ierr,
00015 doublereal *rv4, doublereal *rv5)
00016 {
00017
00018 integer i__1, i__2;
00019 doublereal d__1, d__2, d__3;
00020
00021
00022 static integer i__, j, k, l, p, q, r__, s;
00023 static doublereal u, v;
00024 static integer m1, m2;
00025 static doublereal t1, t2, x0, x1;
00026 static integer ii;
00027 static doublereal xu;
00028 extern doublereal epslon_(doublereal *);
00029 static integer isturm, tag;
00030 static doublereal tst1, tst2;
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
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 --rv5;
00113 --rv4;
00114 --e2;
00115 --e;
00116 --d__;
00117 --ind;
00118 --w;
00119
00120
00121 *ierr = 0;
00122 tag = 0;
00123 t1 = *lb;
00124 t2 = *ub;
00125
00126 i__1 = *n;
00127 for (i__ = 1; i__ <= i__1; ++i__) {
00128 if (i__ == 1) {
00129 goto L20;
00130 }
00131 tst1 = (d__1 = d__[i__], abs(d__1)) + (d__2 = d__[i__ - 1], abs(d__2))
00132 ;
00133 tst2 = tst1 + (d__1 = e[i__], abs(d__1));
00134 if (tst2 > tst1) {
00135 goto L40;
00136 }
00137 L20:
00138 e2[i__] = 0.;
00139 L40:
00140 ;
00141 }
00142
00143
00144 p = 1;
00145 q = *n;
00146 x1 = *ub;
00147 isturm = 1;
00148 goto L320;
00149 L60:
00150 *m = s;
00151 x1 = *lb;
00152 isturm = 2;
00153 goto L320;
00154 L80:
00155 *m -= s;
00156 if (*m > *mm) {
00157 goto L980;
00158 }
00159 q = 0;
00160 r__ = 0;
00161
00162
00163 L100:
00164 if (r__ == *m) {
00165 goto L1001;
00166 }
00167 ++tag;
00168 p = q + 1;
00169 xu = d__[p];
00170 x0 = d__[p];
00171 u = 0.;
00172
00173 i__1 = *n;
00174 for (q = p; q <= i__1; ++q) {
00175 x1 = u;
00176 u = 0.;
00177 v = 0.;
00178 if (q == *n) {
00179 goto L110;
00180 }
00181 u = (d__1 = e[q + 1], abs(d__1));
00182 v = e2[q + 1];
00183 L110:
00184
00185 d__1 = d__[q] - (x1 + u);
00186 xu = min(d__1,xu);
00187
00188 d__1 = d__[q] + (x1 + u);
00189 x0 = max(d__1,x0);
00190 if (v == 0.) {
00191 goto L140;
00192 }
00193
00194 }
00195
00196 L140:
00197
00198 d__2 = abs(xu), d__3 = abs(x0);
00199 d__1 = max(d__2,d__3);
00200 x1 = epslon_(&d__1);
00201 if (*eps1 <= 0.) {
00202 *eps1 = -x1;
00203 }
00204 if (p != q) {
00205 goto L180;
00206 }
00207
00208 if (t1 > d__[p] || d__[p] >= t2) {
00209 goto L940;
00210 }
00211 m1 = p;
00212 m2 = p;
00213 rv5[p] = d__[p];
00214 goto L900;
00215 L180:
00216 x1 *= q - p + 1;
00217
00218 d__1 = t1, d__2 = xu - x1;
00219 *lb = max(d__1,d__2);
00220
00221 d__1 = t2, d__2 = x0 + x1;
00222 *ub = min(d__1,d__2);
00223 x1 = *lb;
00224 isturm = 3;
00225 goto L320;
00226 L200:
00227 m1 = s + 1;
00228 x1 = *ub;
00229 isturm = 4;
00230 goto L320;
00231 L220:
00232 m2 = s;
00233 if (m1 > m2) {
00234 goto L940;
00235 }
00236
00237 x0 = *ub;
00238 isturm = 5;
00239
00240 i__1 = m2;
00241 for (i__ = m1; i__ <= i__1; ++i__) {
00242 rv5[i__] = *ub;
00243 rv4[i__] = *lb;
00244
00245 }
00246
00247
00248
00249
00250 k = m2;
00251 L250:
00252 xu = *lb;
00253
00254 i__1 = k;
00255 for (ii = m1; ii <= i__1; ++ii) {
00256 i__ = m1 + k - ii;
00257 if (xu >= rv4[i__]) {
00258 goto L260;
00259 }
00260 xu = rv4[i__];
00261 goto L280;
00262 L260:
00263 ;
00264 }
00265
00266 L280:
00267 if (x0 > rv5[k]) {
00268 x0 = rv5[k];
00269 }
00270
00271 L300:
00272 x1 = (xu + x0) * .5;
00273 if (x0 - xu <= abs(*eps1)) {
00274 goto L420;
00275 }
00276 tst1 = (abs(xu) + abs(x0)) * 2.;
00277 tst2 = tst1 + (x0 - xu);
00278 if (tst2 == tst1) {
00279 goto L420;
00280 }
00281
00282 L320:
00283 s = p - 1;
00284 u = 1.;
00285
00286 i__1 = q;
00287 for (i__ = p; i__ <= i__1; ++i__) {
00288 if (u != 0.) {
00289 goto L325;
00290 }
00291 v = (d__1 = e[i__], abs(d__1)) / epslon_(&c_b26);
00292 if (e2[i__] == 0.) {
00293 v = 0.;
00294 }
00295 goto L330;
00296 L325:
00297 v = e2[i__] / u;
00298 L330:
00299 u = d__[i__] - x1 - v;
00300 if (u < 0.) {
00301 ++s;
00302 }
00303
00304 }
00305
00306 switch (isturm) {
00307 case 1: goto L60;
00308 case 2: goto L80;
00309 case 3: goto L200;
00310 case 4: goto L220;
00311 case 5: goto L360;
00312 }
00313
00314 L360:
00315 if (s >= k) {
00316 goto L400;
00317 }
00318 xu = x1;
00319 if (s >= m1) {
00320 goto L380;
00321 }
00322 rv4[m1] = x1;
00323 goto L300;
00324 L380:
00325 rv4[s + 1] = x1;
00326 if (rv5[s] > x1) {
00327 rv5[s] = x1;
00328 }
00329 goto L300;
00330 L400:
00331 x0 = x1;
00332 goto L300;
00333
00334 L420:
00335 rv5[k] = x1;
00336 --k;
00337 if (k >= m1) {
00338 goto L250;
00339 }
00340
00341
00342 L900:
00343 s = r__;
00344 r__ = r__ + m2 - m1 + 1;
00345 j = 1;
00346 k = m1;
00347
00348 i__1 = r__;
00349 for (l = 1; l <= i__1; ++l) {
00350 if (j > s) {
00351 goto L910;
00352 }
00353 if (k > m2) {
00354 goto L940;
00355 }
00356 if (rv5[k] >= w[l]) {
00357 goto L915;
00358 }
00359
00360 i__2 = s;
00361 for (ii = j; ii <= i__2; ++ii) {
00362 i__ = l + s - ii;
00363 w[i__ + 1] = w[i__];
00364 ind[i__ + 1] = ind[i__];
00365
00366 }
00367
00368 L910:
00369 w[l] = rv5[k];
00370 ind[l] = tag;
00371 ++k;
00372 goto L920;
00373 L915:
00374 ++j;
00375 L920:
00376 ;
00377 }
00378
00379 L940:
00380 if (q < *n) {
00381 goto L100;
00382 }
00383 goto L1001;
00384
00385
00386 L980:
00387 *ierr = *n * 3 + 1;
00388 L1001:
00389 *lb = t1;
00390 *ub = t2;
00391 return 0;
00392 }
00393