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