00001 /* slaebz.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Subroutine */ int slaebz_(integer *ijob, integer *nitmax, integer *n, 00017 integer *mmax, integer *minp, integer *nbmin, real *abstol, real * 00018 reltol, real *pivmin, real *d__, real *e, real *e2, integer *nval, 00019 real *ab, real *c__, integer *mout, integer *nab, real *work, integer 00020 *iwork, integer *info) 00021 { 00022 /* System generated locals */ 00023 integer nab_dim1, nab_offset, ab_dim1, ab_offset, i__1, i__2, i__3, i__4, 00024 i__5, i__6; 00025 real r__1, r__2, r__3, r__4; 00026 00027 /* Local variables */ 00028 integer j, kf, ji, kl, jp, jit; 00029 real tmp1, tmp2; 00030 integer itmp1, itmp2, kfnew, klnew; 00031 00032 00033 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00034 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00035 /* November 2006 */ 00036 00037 /* .. Scalar Arguments .. */ 00038 /* .. */ 00039 /* .. Array Arguments .. */ 00040 /* .. */ 00041 00042 /* Purpose */ 00043 /* ======= */ 00044 00045 /* SLAEBZ contains the iteration loops which compute and use the */ 00046 /* function N(w), which is the count of eigenvalues of a symmetric */ 00047 /* tridiagonal matrix T less than or equal to its argument w. It */ 00048 /* performs a choice of two types of loops: */ 00049 00050 /* IJOB=1, followed by */ 00051 /* IJOB=2: It takes as input a list of intervals and returns a list of */ 00052 /* sufficiently small intervals whose union contains the same */ 00053 /* eigenvalues as the union of the original intervals. */ 00054 /* The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. */ 00055 /* The output interval (AB(j,1),AB(j,2)] will contain */ 00056 /* eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. */ 00057 00058 /* IJOB=3: It performs a binary search in each input interval */ 00059 /* (AB(j,1),AB(j,2)] for a point w(j) such that */ 00060 /* N(w(j))=NVAL(j), and uses C(j) as the starting point of */ 00061 /* the search. If such a w(j) is found, then on output */ 00062 /* AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output */ 00063 /* (AB(j,1),AB(j,2)] will be a small interval containing the */ 00064 /* point where N(w) jumps through NVAL(j), unless that point */ 00065 /* lies outside the initial interval. */ 00066 00067 /* Note that the intervals are in all cases half-open intervals, */ 00068 /* i.e., of the form (a,b] , which includes b but not a . */ 00069 00070 /* To avoid underflow, the matrix should be scaled so that its largest */ 00071 /* element is no greater than overflow**(1/2) * underflow**(1/4) */ 00072 /* in absolute value. To assure the most accurate computation */ 00073 /* of small eigenvalues, the matrix should be scaled to be */ 00074 /* not much smaller than that, either. */ 00075 00076 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */ 00077 /* Matrix", Report CS41, Computer Science Dept., Stanford */ 00078 /* University, July 21, 1966 */ 00079 00080 /* Note: the arguments are, in general, *not* checked for unreasonable */ 00081 /* values. */ 00082 00083 /* Arguments */ 00084 /* ========= */ 00085 00086 /* IJOB (input) INTEGER */ 00087 /* Specifies what is to be done: */ 00088 /* = 1: Compute NAB for the initial intervals. */ 00089 /* = 2: Perform bisection iteration to find eigenvalues of T. */ 00090 /* = 3: Perform bisection iteration to invert N(w), i.e., */ 00091 /* to find a point which has a specified number of */ 00092 /* eigenvalues of T to its left. */ 00093 /* Other values will cause SLAEBZ to return with INFO=-1. */ 00094 00095 /* NITMAX (input) INTEGER */ 00096 /* The maximum number of "levels" of bisection to be */ 00097 /* performed, i.e., an interval of width W will not be made */ 00098 /* smaller than 2^(-NITMAX) * W. If not all intervals */ 00099 /* have converged after NITMAX iterations, then INFO is set */ 00100 /* to the number of non-converged intervals. */ 00101 00102 /* N (input) INTEGER */ 00103 /* The dimension n of the tridiagonal matrix T. It must be at */ 00104 /* least 1. */ 00105 00106 /* MMAX (input) INTEGER */ 00107 /* The maximum number of intervals. If more than MMAX intervals */ 00108 /* are generated, then SLAEBZ will quit with INFO=MMAX+1. */ 00109 00110 /* MINP (input) INTEGER */ 00111 /* The initial number of intervals. It may not be greater than */ 00112 /* MMAX. */ 00113 00114 /* NBMIN (input) INTEGER */ 00115 /* The smallest number of intervals that should be processed */ 00116 /* using a vector loop. If zero, then only the scalar loop */ 00117 /* will be used. */ 00118 00119 /* ABSTOL (input) REAL */ 00120 /* The minimum (absolute) width of an interval. When an */ 00121 /* interval is narrower than ABSTOL, or than RELTOL times the */ 00122 /* larger (in magnitude) endpoint, then it is considered to be */ 00123 /* sufficiently small, i.e., converged. This must be at least */ 00124 /* zero. */ 00125 00126 /* RELTOL (input) REAL */ 00127 /* The minimum relative width of an interval. When an interval */ 00128 /* is narrower than ABSTOL, or than RELTOL times the larger (in */ 00129 /* magnitude) endpoint, then it is considered to be */ 00130 /* sufficiently small, i.e., converged. Note: this should */ 00131 /* always be at least radix*machine epsilon. */ 00132 00133 /* PIVMIN (input) REAL */ 00134 /* The minimum absolute value of a "pivot" in the Sturm */ 00135 /* sequence loop. This *must* be at least max |e(j)**2| * */ 00136 /* safe_min and at least safe_min, where safe_min is at least */ 00137 /* the smallest number that can divide one without overflow. */ 00138 00139 /* D (input) REAL array, dimension (N) */ 00140 /* The diagonal elements of the tridiagonal matrix T. */ 00141 00142 /* E (input) REAL array, dimension (N) */ 00143 /* The offdiagonal elements of the tridiagonal matrix T in */ 00144 /* positions 1 through N-1. E(N) is arbitrary. */ 00145 00146 /* E2 (input) REAL array, dimension (N) */ 00147 /* The squares of the offdiagonal elements of the tridiagonal */ 00148 /* matrix T. E2(N) is ignored. */ 00149 00150 /* NVAL (input/output) INTEGER array, dimension (MINP) */ 00151 /* If IJOB=1 or 2, not referenced. */ 00152 /* If IJOB=3, the desired values of N(w). The elements of NVAL */ 00153 /* will be reordered to correspond with the intervals in AB. */ 00154 /* Thus, NVAL(j) on output will not, in general be the same as */ 00155 /* NVAL(j) on input, but it will correspond with the interval */ 00156 /* (AB(j,1),AB(j,2)] on output. */ 00157 00158 /* AB (input/output) REAL array, dimension (MMAX,2) */ 00159 /* The endpoints of the intervals. AB(j,1) is a(j), the left */ 00160 /* endpoint of the j-th interval, and AB(j,2) is b(j), the */ 00161 /* right endpoint of the j-th interval. The input intervals */ 00162 /* will, in general, be modified, split, and reordered by the */ 00163 /* calculation. */ 00164 00165 /* C (input/output) REAL array, dimension (MMAX) */ 00166 /* If IJOB=1, ignored. */ 00167 /* If IJOB=2, workspace. */ 00168 /* If IJOB=3, then on input C(j) should be initialized to the */ 00169 /* first search point in the binary search. */ 00170 00171 /* MOUT (output) INTEGER */ 00172 /* If IJOB=1, the number of eigenvalues in the intervals. */ 00173 /* If IJOB=2 or 3, the number of intervals output. */ 00174 /* If IJOB=3, MOUT will equal MINP. */ 00175 00176 /* NAB (input/output) INTEGER array, dimension (MMAX,2) */ 00177 /* If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). */ 00178 /* If IJOB=2, then on input, NAB(i,j) should be set. It must */ 00179 /* satisfy the condition: */ 00180 /* N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), */ 00181 /* which means that in interval i only eigenvalues */ 00182 /* NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, */ 00183 /* NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with */ 00184 /* IJOB=1. */ 00185 /* On output, NAB(i,j) will contain */ 00186 /* max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of */ 00187 /* the input interval that the output interval */ 00188 /* (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the */ 00189 /* the input values of NAB(k,1) and NAB(k,2). */ 00190 /* If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), */ 00191 /* unless N(w) > NVAL(i) for all search points w , in which */ 00192 /* case NAB(i,1) will not be modified, i.e., the output */ 00193 /* value will be the same as the input value (modulo */ 00194 /* reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) */ 00195 /* for all search points w , in which case NAB(i,2) will */ 00196 /* not be modified. Normally, NAB should be set to some */ 00197 /* distinctive value(s) before SLAEBZ is called. */ 00198 00199 /* WORK (workspace) REAL array, dimension (MMAX) */ 00200 /* Workspace. */ 00201 00202 /* IWORK (workspace) INTEGER array, dimension (MMAX) */ 00203 /* Workspace. */ 00204 00205 /* INFO (output) INTEGER */ 00206 /* = 0: All intervals converged. */ 00207 /* = 1--MMAX: The last INFO intervals did not converge. */ 00208 /* = MMAX+1: More than MMAX intervals were generated. */ 00209 00210 /* Further Details */ 00211 /* =============== */ 00212 00213 /* This routine is intended to be called only by other LAPACK */ 00214 /* routines, thus the interface is less user-friendly. It is intended */ 00215 /* for two purposes: */ 00216 00217 /* (a) finding eigenvalues. In this case, SLAEBZ should have one or */ 00218 /* more initial intervals set up in AB, and SLAEBZ should be called */ 00219 /* with IJOB=1. This sets up NAB, and also counts the eigenvalues. */ 00220 /* Intervals with no eigenvalues would usually be thrown out at */ 00221 /* this point. Also, if not all the eigenvalues in an interval i */ 00222 /* are desired, NAB(i,1) can be increased or NAB(i,2) decreased. */ 00223 /* For example, set NAB(i,1)=NAB(i,2)-1 to get the largest */ 00224 /* eigenvalue. SLAEBZ is then called with IJOB=2 and MMAX */ 00225 /* no smaller than the value of MOUT returned by the call with */ 00226 /* IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 */ 00227 /* through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the */ 00228 /* tolerance specified by ABSTOL and RELTOL. */ 00229 00230 /* (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). */ 00231 /* In this case, start with a Gershgorin interval (a,b). Set up */ 00232 /* AB to contain 2 search intervals, both initially (a,b). One */ 00233 /* NVAL element should contain f-1 and the other should contain l */ 00234 /* , while C should contain a and b, resp. NAB(i,1) should be -1 */ 00235 /* and NAB(i,2) should be N+1, to flag an error if the desired */ 00236 /* interval does not lie in (a,b). SLAEBZ is then called with */ 00237 /* IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- */ 00238 /* j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while */ 00239 /* if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r */ 00240 /* >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and */ 00241 /* N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and */ 00242 /* w(l-r)=...=w(l+k) are handled similarly. */ 00243 00244 /* ===================================================================== */ 00245 00246 /* .. Parameters .. */ 00247 /* .. */ 00248 /* .. Local Scalars .. */ 00249 /* .. */ 00250 /* .. Intrinsic Functions .. */ 00251 /* .. */ 00252 /* .. Executable Statements .. */ 00253 00254 /* Check for Errors */ 00255 00256 /* Parameter adjustments */ 00257 nab_dim1 = *mmax; 00258 nab_offset = 1 + nab_dim1; 00259 nab -= nab_offset; 00260 ab_dim1 = *mmax; 00261 ab_offset = 1 + ab_dim1; 00262 ab -= ab_offset; 00263 --d__; 00264 --e; 00265 --e2; 00266 --nval; 00267 --c__; 00268 --work; 00269 --iwork; 00270 00271 /* Function Body */ 00272 *info = 0; 00273 if (*ijob < 1 || *ijob > 3) { 00274 *info = -1; 00275 return 0; 00276 } 00277 00278 /* Initialize NAB */ 00279 00280 if (*ijob == 1) { 00281 00282 /* Compute the number of eigenvalues in the initial intervals. */ 00283 00284 *mout = 0; 00285 /* DIR$ NOVECTOR */ 00286 i__1 = *minp; 00287 for (ji = 1; ji <= i__1; ++ji) { 00288 for (jp = 1; jp <= 2; ++jp) { 00289 tmp1 = d__[1] - ab[ji + jp * ab_dim1]; 00290 if (dabs(tmp1) < *pivmin) { 00291 tmp1 = -(*pivmin); 00292 } 00293 nab[ji + jp * nab_dim1] = 0; 00294 if (tmp1 <= 0.f) { 00295 nab[ji + jp * nab_dim1] = 1; 00296 } 00297 00298 i__2 = *n; 00299 for (j = 2; j <= i__2; ++j) { 00300 tmp1 = d__[j] - e2[j - 1] / tmp1 - ab[ji + jp * ab_dim1]; 00301 if (dabs(tmp1) < *pivmin) { 00302 tmp1 = -(*pivmin); 00303 } 00304 if (tmp1 <= 0.f) { 00305 ++nab[ji + jp * nab_dim1]; 00306 } 00307 /* L10: */ 00308 } 00309 /* L20: */ 00310 } 00311 *mout = *mout + nab[ji + (nab_dim1 << 1)] - nab[ji + nab_dim1]; 00312 /* L30: */ 00313 } 00314 return 0; 00315 } 00316 00317 /* Initialize for loop */ 00318 00319 /* KF and KL have the following meaning: */ 00320 /* Intervals 1,...,KF-1 have converged. */ 00321 /* Intervals KF,...,KL still need to be refined. */ 00322 00323 kf = 1; 00324 kl = *minp; 00325 00326 /* If IJOB=2, initialize C. */ 00327 /* If IJOB=3, use the user-supplied starting point. */ 00328 00329 if (*ijob == 2) { 00330 i__1 = *minp; 00331 for (ji = 1; ji <= i__1; ++ji) { 00332 c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5f; 00333 /* L40: */ 00334 } 00335 } 00336 00337 /* Iteration loop */ 00338 00339 i__1 = *nitmax; 00340 for (jit = 1; jit <= i__1; ++jit) { 00341 00342 /* Loop over intervals */ 00343 00344 if (kl - kf + 1 >= *nbmin && *nbmin > 0) { 00345 00346 /* Begin of Parallel Version of the loop */ 00347 00348 i__2 = kl; 00349 for (ji = kf; ji <= i__2; ++ji) { 00350 00351 /* Compute N(c), the number of eigenvalues less than c */ 00352 00353 work[ji] = d__[1] - c__[ji]; 00354 iwork[ji] = 0; 00355 if (work[ji] <= *pivmin) { 00356 iwork[ji] = 1; 00357 /* Computing MIN */ 00358 r__1 = work[ji], r__2 = -(*pivmin); 00359 work[ji] = dmin(r__1,r__2); 00360 } 00361 00362 i__3 = *n; 00363 for (j = 2; j <= i__3; ++j) { 00364 work[ji] = d__[j] - e2[j - 1] / work[ji] - c__[ji]; 00365 if (work[ji] <= *pivmin) { 00366 ++iwork[ji]; 00367 /* Computing MIN */ 00368 r__1 = work[ji], r__2 = -(*pivmin); 00369 work[ji] = dmin(r__1,r__2); 00370 } 00371 /* L50: */ 00372 } 00373 /* L60: */ 00374 } 00375 00376 if (*ijob <= 2) { 00377 00378 /* IJOB=2: Choose all intervals containing eigenvalues. */ 00379 00380 klnew = kl; 00381 i__2 = kl; 00382 for (ji = kf; ji <= i__2; ++ji) { 00383 00384 /* Insure that N(w) is monotone */ 00385 00386 /* Computing MIN */ 00387 /* Computing MAX */ 00388 i__5 = nab[ji + nab_dim1], i__6 = iwork[ji]; 00389 i__3 = nab[ji + (nab_dim1 << 1)], i__4 = max(i__5,i__6); 00390 iwork[ji] = min(i__3,i__4); 00391 00392 /* Update the Queue -- add intervals if both halves */ 00393 /* contain eigenvalues. */ 00394 00395 if (iwork[ji] == nab[ji + (nab_dim1 << 1)]) { 00396 00397 /* No eigenvalue in the upper interval: */ 00398 /* just use the lower interval. */ 00399 00400 ab[ji + (ab_dim1 << 1)] = c__[ji]; 00401 00402 } else if (iwork[ji] == nab[ji + nab_dim1]) { 00403 00404 /* No eigenvalue in the lower interval: */ 00405 /* just use the upper interval. */ 00406 00407 ab[ji + ab_dim1] = c__[ji]; 00408 } else { 00409 ++klnew; 00410 if (klnew <= *mmax) { 00411 00412 /* Eigenvalue in both intervals -- add upper to */ 00413 /* queue. */ 00414 00415 ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 << 00416 1)]; 00417 nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1 00418 << 1)]; 00419 ab[klnew + ab_dim1] = c__[ji]; 00420 nab[klnew + nab_dim1] = iwork[ji]; 00421 ab[ji + (ab_dim1 << 1)] = c__[ji]; 00422 nab[ji + (nab_dim1 << 1)] = iwork[ji]; 00423 } else { 00424 *info = *mmax + 1; 00425 } 00426 } 00427 /* L70: */ 00428 } 00429 if (*info != 0) { 00430 return 0; 00431 } 00432 kl = klnew; 00433 } else { 00434 00435 /* IJOB=3: Binary search. Keep only the interval containing */ 00436 /* w s.t. N(w) = NVAL */ 00437 00438 i__2 = kl; 00439 for (ji = kf; ji <= i__2; ++ji) { 00440 if (iwork[ji] <= nval[ji]) { 00441 ab[ji + ab_dim1] = c__[ji]; 00442 nab[ji + nab_dim1] = iwork[ji]; 00443 } 00444 if (iwork[ji] >= nval[ji]) { 00445 ab[ji + (ab_dim1 << 1)] = c__[ji]; 00446 nab[ji + (nab_dim1 << 1)] = iwork[ji]; 00447 } 00448 /* L80: */ 00449 } 00450 } 00451 00452 } else { 00453 00454 /* End of Parallel Version of the loop */ 00455 00456 /* Begin of Serial Version of the loop */ 00457 00458 klnew = kl; 00459 i__2 = kl; 00460 for (ji = kf; ji <= i__2; ++ji) { 00461 00462 /* Compute N(w), the number of eigenvalues less than w */ 00463 00464 tmp1 = c__[ji]; 00465 tmp2 = d__[1] - tmp1; 00466 itmp1 = 0; 00467 if (tmp2 <= *pivmin) { 00468 itmp1 = 1; 00469 /* Computing MIN */ 00470 r__1 = tmp2, r__2 = -(*pivmin); 00471 tmp2 = dmin(r__1,r__2); 00472 } 00473 00474 /* A series of compiler directives to defeat vectorization */ 00475 /* for the next loop */ 00476 00477 /* $PL$ CMCHAR=' ' */ 00478 /* DIR$ NEXTSCALAR */ 00479 /* $DIR SCALAR */ 00480 /* DIR$ NEXT SCALAR */ 00481 /* VD$L NOVECTOR */ 00482 /* DEC$ NOVECTOR */ 00483 /* VD$ NOVECTOR */ 00484 /* VDIR NOVECTOR */ 00485 /* VOCL LOOP,SCALAR */ 00486 /* IBM PREFER SCALAR */ 00487 /* $PL$ CMCHAR='*' */ 00488 00489 i__3 = *n; 00490 for (j = 2; j <= i__3; ++j) { 00491 tmp2 = d__[j] - e2[j - 1] / tmp2 - tmp1; 00492 if (tmp2 <= *pivmin) { 00493 ++itmp1; 00494 /* Computing MIN */ 00495 r__1 = tmp2, r__2 = -(*pivmin); 00496 tmp2 = dmin(r__1,r__2); 00497 } 00498 /* L90: */ 00499 } 00500 00501 if (*ijob <= 2) { 00502 00503 /* IJOB=2: Choose all intervals containing eigenvalues. */ 00504 00505 /* Insure that N(w) is monotone */ 00506 00507 /* Computing MIN */ 00508 /* Computing MAX */ 00509 i__5 = nab[ji + nab_dim1]; 00510 i__3 = nab[ji + (nab_dim1 << 1)], i__4 = max(i__5,itmp1); 00511 itmp1 = min(i__3,i__4); 00512 00513 /* Update the Queue -- add intervals if both halves */ 00514 /* contain eigenvalues. */ 00515 00516 if (itmp1 == nab[ji + (nab_dim1 << 1)]) { 00517 00518 /* No eigenvalue in the upper interval: */ 00519 /* just use the lower interval. */ 00520 00521 ab[ji + (ab_dim1 << 1)] = tmp1; 00522 00523 } else if (itmp1 == nab[ji + nab_dim1]) { 00524 00525 /* No eigenvalue in the lower interval: */ 00526 /* just use the upper interval. */ 00527 00528 ab[ji + ab_dim1] = tmp1; 00529 } else if (klnew < *mmax) { 00530 00531 /* Eigenvalue in both intervals -- add upper to queue. */ 00532 00533 ++klnew; 00534 ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 << 1)]; 00535 nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1 << 00536 1)]; 00537 ab[klnew + ab_dim1] = tmp1; 00538 nab[klnew + nab_dim1] = itmp1; 00539 ab[ji + (ab_dim1 << 1)] = tmp1; 00540 nab[ji + (nab_dim1 << 1)] = itmp1; 00541 } else { 00542 *info = *mmax + 1; 00543 return 0; 00544 } 00545 } else { 00546 00547 /* IJOB=3: Binary search. Keep only the interval */ 00548 /* containing w s.t. N(w) = NVAL */ 00549 00550 if (itmp1 <= nval[ji]) { 00551 ab[ji + ab_dim1] = tmp1; 00552 nab[ji + nab_dim1] = itmp1; 00553 } 00554 if (itmp1 >= nval[ji]) { 00555 ab[ji + (ab_dim1 << 1)] = tmp1; 00556 nab[ji + (nab_dim1 << 1)] = itmp1; 00557 } 00558 } 00559 /* L100: */ 00560 } 00561 kl = klnew; 00562 00563 /* End of Serial Version of the loop */ 00564 00565 } 00566 00567 /* Check for convergence */ 00568 00569 kfnew = kf; 00570 i__2 = kl; 00571 for (ji = kf; ji <= i__2; ++ji) { 00572 tmp1 = (r__1 = ab[ji + (ab_dim1 << 1)] - ab[ji + ab_dim1], dabs( 00573 r__1)); 00574 /* Computing MAX */ 00575 r__3 = (r__1 = ab[ji + (ab_dim1 << 1)], dabs(r__1)), r__4 = (r__2 00576 = ab[ji + ab_dim1], dabs(r__2)); 00577 tmp2 = dmax(r__3,r__4); 00578 /* Computing MAX */ 00579 r__1 = max(*abstol,*pivmin), r__2 = *reltol * tmp2; 00580 if (tmp1 < dmax(r__1,r__2) || nab[ji + nab_dim1] >= nab[ji + ( 00581 nab_dim1 << 1)]) { 00582 00583 /* Converged -- Swap with position KFNEW, */ 00584 /* then increment KFNEW */ 00585 00586 if (ji > kfnew) { 00587 tmp1 = ab[ji + ab_dim1]; 00588 tmp2 = ab[ji + (ab_dim1 << 1)]; 00589 itmp1 = nab[ji + nab_dim1]; 00590 itmp2 = nab[ji + (nab_dim1 << 1)]; 00591 ab[ji + ab_dim1] = ab[kfnew + ab_dim1]; 00592 ab[ji + (ab_dim1 << 1)] = ab[kfnew + (ab_dim1 << 1)]; 00593 nab[ji + nab_dim1] = nab[kfnew + nab_dim1]; 00594 nab[ji + (nab_dim1 << 1)] = nab[kfnew + (nab_dim1 << 1)]; 00595 ab[kfnew + ab_dim1] = tmp1; 00596 ab[kfnew + (ab_dim1 << 1)] = tmp2; 00597 nab[kfnew + nab_dim1] = itmp1; 00598 nab[kfnew + (nab_dim1 << 1)] = itmp2; 00599 if (*ijob == 3) { 00600 itmp1 = nval[ji]; 00601 nval[ji] = nval[kfnew]; 00602 nval[kfnew] = itmp1; 00603 } 00604 } 00605 ++kfnew; 00606 } 00607 /* L110: */ 00608 } 00609 kf = kfnew; 00610 00611 /* Choose Midpoints */ 00612 00613 i__2 = kl; 00614 for (ji = kf; ji <= i__2; ++ji) { 00615 c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5f; 00616 /* L120: */ 00617 } 00618 00619 /* If no more intervals to refine, quit. */ 00620 00621 if (kf > kl) { 00622 goto L140; 00623 } 00624 /* L130: */ 00625 } 00626 00627 /* Converged */ 00628 00629 L140: 00630 /* Computing MAX */ 00631 i__1 = kl + 1 - kf; 00632 *info = max(i__1,0); 00633 *mout = kl; 00634 00635 return 0; 00636 00637 /* End of SLAEBZ */ 00638 00639 } /* slaebz_ */