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