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