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