1*> \brief \b DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLAEBZ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaebz.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaebz.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaebz.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, 22* RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, 23* NAB, WORK, IWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX 27* DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL 28* .. 29* .. Array Arguments .. 30* INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 31* DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 32* $ WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> DLAEBZ contains the iteration loops which compute and use the 42*> function N(w), which is the count of eigenvalues of a symmetric 43*> tridiagonal matrix T less than or equal to its argument w. It 44*> performs a choice of two types of loops: 45*> 46*> IJOB=1, followed by 47*> IJOB=2: It takes as input a list of intervals and returns a list of 48*> sufficiently small intervals whose union contains the same 49*> eigenvalues as the union of the original intervals. 50*> The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. 51*> The output interval (AB(j,1),AB(j,2)] will contain 52*> eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. 53*> 54*> IJOB=3: It performs a binary search in each input interval 55*> (AB(j,1),AB(j,2)] for a point w(j) such that 56*> N(w(j))=NVAL(j), and uses C(j) as the starting point of 57*> the search. If such a w(j) is found, then on output 58*> AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output 59*> (AB(j,1),AB(j,2)] will be a small interval containing the 60*> point where N(w) jumps through NVAL(j), unless that point 61*> lies outside the initial interval. 62*> 63*> Note that the intervals are in all cases half-open intervals, 64*> i.e., of the form (a,b] , which includes b but not a . 65*> 66*> To avoid underflow, the matrix should be scaled so that its largest 67*> element is no greater than overflow**(1/2) * underflow**(1/4) 68*> in absolute value. To assure the most accurate computation 69*> of small eigenvalues, the matrix should be scaled to be 70*> not much smaller than that, either. 71*> 72*> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 73*> Matrix", Report CS41, Computer Science Dept., Stanford 74*> University, July 21, 1966 75*> 76*> Note: the arguments are, in general, *not* checked for unreasonable 77*> values. 78*> \endverbatim 79* 80* Arguments: 81* ========== 82* 83*> \param[in] IJOB 84*> \verbatim 85*> IJOB is INTEGER 86*> Specifies what is to be done: 87*> = 1: Compute NAB for the initial intervals. 88*> = 2: Perform bisection iteration to find eigenvalues of T. 89*> = 3: Perform bisection iteration to invert N(w), i.e., 90*> to find a point which has a specified number of 91*> eigenvalues of T to its left. 92*> Other values will cause DLAEBZ to return with INFO=-1. 93*> \endverbatim 94*> 95*> \param[in] NITMAX 96*> \verbatim 97*> NITMAX is INTEGER 98*> The maximum number of "levels" of bisection to be 99*> performed, i.e., an interval of width W will not be made 100*> smaller than 2^(-NITMAX) * W. If not all intervals 101*> have converged after NITMAX iterations, then INFO is set 102*> to the number of non-converged intervals. 103*> \endverbatim 104*> 105*> \param[in] N 106*> \verbatim 107*> N is INTEGER 108*> The dimension n of the tridiagonal matrix T. It must be at 109*> least 1. 110*> \endverbatim 111*> 112*> \param[in] MMAX 113*> \verbatim 114*> MMAX is INTEGER 115*> The maximum number of intervals. If more than MMAX intervals 116*> are generated, then DLAEBZ will quit with INFO=MMAX+1. 117*> \endverbatim 118*> 119*> \param[in] MINP 120*> \verbatim 121*> MINP is INTEGER 122*> The initial number of intervals. It may not be greater than 123*> MMAX. 124*> \endverbatim 125*> 126*> \param[in] NBMIN 127*> \verbatim 128*> NBMIN is INTEGER 129*> The smallest number of intervals that should be processed 130*> using a vector loop. If zero, then only the scalar loop 131*> will be used. 132*> \endverbatim 133*> 134*> \param[in] ABSTOL 135*> \verbatim 136*> ABSTOL is DOUBLE PRECISION 137*> The minimum (absolute) width of an interval. When an 138*> interval is narrower than ABSTOL, or than RELTOL times the 139*> larger (in magnitude) endpoint, then it is considered to be 140*> sufficiently small, i.e., converged. This must be at least 141*> zero. 142*> \endverbatim 143*> 144*> \param[in] RELTOL 145*> \verbatim 146*> RELTOL is DOUBLE PRECISION 147*> The minimum relative width of an interval. When an interval 148*> is narrower than ABSTOL, or than RELTOL times the larger (in 149*> magnitude) endpoint, then it is considered to be 150*> sufficiently small, i.e., converged. Note: this should 151*> always be at least radix*machine epsilon. 152*> \endverbatim 153*> 154*> \param[in] PIVMIN 155*> \verbatim 156*> PIVMIN is DOUBLE PRECISION 157*> The minimum absolute value of a "pivot" in the Sturm 158*> sequence loop. 159*> This must be at least max |e(j)**2|*safe_min and at 160*> least safe_min, where safe_min is at least 161*> the smallest number that can divide one without overflow. 162*> \endverbatim 163*> 164*> \param[in] D 165*> \verbatim 166*> D is DOUBLE PRECISION array, dimension (N) 167*> The diagonal elements of the tridiagonal matrix T. 168*> \endverbatim 169*> 170*> \param[in] E 171*> \verbatim 172*> E is DOUBLE PRECISION array, dimension (N) 173*> The offdiagonal elements of the tridiagonal matrix T in 174*> positions 1 through N-1. E(N) is arbitrary. 175*> \endverbatim 176*> 177*> \param[in] E2 178*> \verbatim 179*> E2 is DOUBLE PRECISION array, dimension (N) 180*> The squares of the offdiagonal elements of the tridiagonal 181*> matrix T. E2(N) is ignored. 182*> \endverbatim 183*> 184*> \param[in,out] NVAL 185*> \verbatim 186*> NVAL is INTEGER array, dimension (MINP) 187*> If IJOB=1 or 2, not referenced. 188*> If IJOB=3, the desired values of N(w). The elements of NVAL 189*> will be reordered to correspond with the intervals in AB. 190*> Thus, NVAL(j) on output will not, in general be the same as 191*> NVAL(j) on input, but it will correspond with the interval 192*> (AB(j,1),AB(j,2)] on output. 193*> \endverbatim 194*> 195*> \param[in,out] AB 196*> \verbatim 197*> AB is DOUBLE PRECISION array, dimension (MMAX,2) 198*> The endpoints of the intervals. AB(j,1) is a(j), the left 199*> endpoint of the j-th interval, and AB(j,2) is b(j), the 200*> right endpoint of the j-th interval. The input intervals 201*> will, in general, be modified, split, and reordered by the 202*> calculation. 203*> \endverbatim 204*> 205*> \param[in,out] C 206*> \verbatim 207*> C is DOUBLE PRECISION array, dimension (MMAX) 208*> If IJOB=1, ignored. 209*> If IJOB=2, workspace. 210*> If IJOB=3, then on input C(j) should be initialized to the 211*> first search point in the binary search. 212*> \endverbatim 213*> 214*> \param[out] MOUT 215*> \verbatim 216*> MOUT is INTEGER 217*> If IJOB=1, the number of eigenvalues in the intervals. 218*> If IJOB=2 or 3, the number of intervals output. 219*> If IJOB=3, MOUT will equal MINP. 220*> \endverbatim 221*> 222*> \param[in,out] NAB 223*> \verbatim 224*> NAB is INTEGER array, dimension (MMAX,2) 225*> If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). 226*> If IJOB=2, then on input, NAB(i,j) should be set. It must 227*> satisfy the condition: 228*> N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), 229*> which means that in interval i only eigenvalues 230*> NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, 231*> NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with 232*> IJOB=1. 233*> On output, NAB(i,j) will contain 234*> max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of 235*> the input interval that the output interval 236*> (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the 237*> the input values of NAB(k,1) and NAB(k,2). 238*> If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), 239*> unless N(w) > NVAL(i) for all search points w , in which 240*> case NAB(i,1) will not be modified, i.e., the output 241*> value will be the same as the input value (modulo 242*> reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) 243*> for all search points w , in which case NAB(i,2) will 244*> not be modified. Normally, NAB should be set to some 245*> distinctive value(s) before DLAEBZ is called. 246*> \endverbatim 247*> 248*> \param[out] WORK 249*> \verbatim 250*> WORK is DOUBLE PRECISION array, dimension (MMAX) 251*> Workspace. 252*> \endverbatim 253*> 254*> \param[out] IWORK 255*> \verbatim 256*> IWORK is INTEGER array, dimension (MMAX) 257*> Workspace. 258*> \endverbatim 259*> 260*> \param[out] INFO 261*> \verbatim 262*> INFO is INTEGER 263*> = 0: All intervals converged. 264*> = 1--MMAX: The last INFO intervals did not converge. 265*> = MMAX+1: More than MMAX intervals were generated. 266*> \endverbatim 267* 268* Authors: 269* ======== 270* 271*> \author Univ. of Tennessee 272*> \author Univ. of California Berkeley 273*> \author Univ. of Colorado Denver 274*> \author NAG Ltd. 275* 276*> \ingroup OTHERauxiliary 277* 278*> \par Further Details: 279* ===================== 280*> 281*> \verbatim 282*> 283*> This routine is intended to be called only by other LAPACK 284*> routines, thus the interface is less user-friendly. It is intended 285*> for two purposes: 286*> 287*> (a) finding eigenvalues. In this case, DLAEBZ should have one or 288*> more initial intervals set up in AB, and DLAEBZ should be called 289*> with IJOB=1. This sets up NAB, and also counts the eigenvalues. 290*> Intervals with no eigenvalues would usually be thrown out at 291*> this point. Also, if not all the eigenvalues in an interval i 292*> are desired, NAB(i,1) can be increased or NAB(i,2) decreased. 293*> For example, set NAB(i,1)=NAB(i,2)-1 to get the largest 294*> eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX 295*> no smaller than the value of MOUT returned by the call with 296*> IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 297*> through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the 298*> tolerance specified by ABSTOL and RELTOL. 299*> 300*> (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). 301*> In this case, start with a Gershgorin interval (a,b). Set up 302*> AB to contain 2 search intervals, both initially (a,b). One 303*> NVAL element should contain f-1 and the other should contain l 304*> , while C should contain a and b, resp. NAB(i,1) should be -1 305*> and NAB(i,2) should be N+1, to flag an error if the desired 306*> interval does not lie in (a,b). DLAEBZ is then called with 307*> IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- 308*> j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while 309*> if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r 310*> >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and 311*> N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and 312*> w(l-r)=...=w(l+k) are handled similarly. 313*> \endverbatim 314*> 315* ===================================================================== 316 SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, 317 $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, 318 $ NAB, WORK, IWORK, INFO ) 319* 320* -- LAPACK auxiliary routine -- 321* -- LAPACK is a software package provided by Univ. of Tennessee, -- 322* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 323* 324* .. Scalar Arguments .. 325 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX 326 DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL 327* .. 328* .. Array Arguments .. 329 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 330 DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 331 $ WORK( * ) 332* .. 333* 334* ===================================================================== 335* 336* .. Parameters .. 337 DOUBLE PRECISION ZERO, TWO, HALF 338 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, 339 $ HALF = 1.0D0 / TWO ) 340* .. 341* .. Local Scalars .. 342 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, 343 $ KLNEW 344 DOUBLE PRECISION TMP1, TMP2 345* .. 346* .. Intrinsic Functions .. 347 INTRINSIC ABS, MAX, MIN 348* .. 349* .. Executable Statements .. 350* 351* Check for Errors 352* 353 INFO = 0 354 IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN 355 INFO = -1 356 RETURN 357 END IF 358* 359* Initialize NAB 360* 361 IF( IJOB.EQ.1 ) THEN 362* 363* Compute the number of eigenvalues in the initial intervals. 364* 365 MOUT = 0 366 DO 30 JI = 1, MINP 367 DO 20 JP = 1, 2 368 TMP1 = D( 1 ) - AB( JI, JP ) 369 IF( ABS( TMP1 ).LT.PIVMIN ) 370 $ TMP1 = -PIVMIN 371 NAB( JI, JP ) = 0 372 IF( TMP1.LE.ZERO ) 373 $ NAB( JI, JP ) = 1 374* 375 DO 10 J = 2, N 376 TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) 377 IF( ABS( TMP1 ).LT.PIVMIN ) 378 $ TMP1 = -PIVMIN 379 IF( TMP1.LE.ZERO ) 380 $ NAB( JI, JP ) = NAB( JI, JP ) + 1 381 10 CONTINUE 382 20 CONTINUE 383 MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) 384 30 CONTINUE 385 RETURN 386 END IF 387* 388* Initialize for loop 389* 390* KF and KL have the following meaning: 391* Intervals 1,...,KF-1 have converged. 392* Intervals KF,...,KL still need to be refined. 393* 394 KF = 1 395 KL = MINP 396* 397* If IJOB=2, initialize C. 398* If IJOB=3, use the user-supplied starting point. 399* 400 IF( IJOB.EQ.2 ) THEN 401 DO 40 JI = 1, MINP 402 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 403 40 CONTINUE 404 END IF 405* 406* Iteration loop 407* 408 DO 130 JIT = 1, NITMAX 409* 410* Loop over intervals 411* 412 IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN 413* 414* Begin of Parallel Version of the loop 415* 416 DO 60 JI = KF, KL 417* 418* Compute N(c), the number of eigenvalues less than c 419* 420 WORK( JI ) = D( 1 ) - C( JI ) 421 IWORK( JI ) = 0 422 IF( WORK( JI ).LE.PIVMIN ) THEN 423 IWORK( JI ) = 1 424 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 425 END IF 426* 427 DO 50 J = 2, N 428 WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) 429 IF( WORK( JI ).LE.PIVMIN ) THEN 430 IWORK( JI ) = IWORK( JI ) + 1 431 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 432 END IF 433 50 CONTINUE 434 60 CONTINUE 435* 436 IF( IJOB.LE.2 ) THEN 437* 438* IJOB=2: Choose all intervals containing eigenvalues. 439* 440 KLNEW = KL 441 DO 70 JI = KF, KL 442* 443* Insure that N(w) is monotone 444* 445 IWORK( JI ) = MIN( NAB( JI, 2 ), 446 $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) 447* 448* Update the Queue -- add intervals if both halves 449* contain eigenvalues. 450* 451 IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN 452* 453* No eigenvalue in the upper interval: 454* just use the lower interval. 455* 456 AB( JI, 2 ) = C( JI ) 457* 458 ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN 459* 460* No eigenvalue in the lower interval: 461* just use the upper interval. 462* 463 AB( JI, 1 ) = C( JI ) 464 ELSE 465 KLNEW = KLNEW + 1 466 IF( KLNEW.LE.MMAX ) THEN 467* 468* Eigenvalue in both intervals -- add upper to 469* queue. 470* 471 AB( KLNEW, 2 ) = AB( JI, 2 ) 472 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 473 AB( KLNEW, 1 ) = C( JI ) 474 NAB( KLNEW, 1 ) = IWORK( JI ) 475 AB( JI, 2 ) = C( JI ) 476 NAB( JI, 2 ) = IWORK( JI ) 477 ELSE 478 INFO = MMAX + 1 479 END IF 480 END IF 481 70 CONTINUE 482 IF( INFO.NE.0 ) 483 $ RETURN 484 KL = KLNEW 485 ELSE 486* 487* IJOB=3: Binary search. Keep only the interval containing 488* w s.t. N(w) = NVAL 489* 490 DO 80 JI = KF, KL 491 IF( IWORK( JI ).LE.NVAL( JI ) ) THEN 492 AB( JI, 1 ) = C( JI ) 493 NAB( JI, 1 ) = IWORK( JI ) 494 END IF 495 IF( IWORK( JI ).GE.NVAL( JI ) ) THEN 496 AB( JI, 2 ) = C( JI ) 497 NAB( JI, 2 ) = IWORK( JI ) 498 END IF 499 80 CONTINUE 500 END IF 501* 502 ELSE 503* 504* End of Parallel Version of the loop 505* 506* Begin of Serial Version of the loop 507* 508 KLNEW = KL 509 DO 100 JI = KF, KL 510* 511* Compute N(w), the number of eigenvalues less than w 512* 513 TMP1 = C( JI ) 514 TMP2 = D( 1 ) - TMP1 515 ITMP1 = 0 516 IF( TMP2.LE.PIVMIN ) THEN 517 ITMP1 = 1 518 TMP2 = MIN( TMP2, -PIVMIN ) 519 END IF 520* 521 DO 90 J = 2, N 522 TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 523 IF( TMP2.LE.PIVMIN ) THEN 524 ITMP1 = ITMP1 + 1 525 TMP2 = MIN( TMP2, -PIVMIN ) 526 END IF 527 90 CONTINUE 528* 529 IF( IJOB.LE.2 ) THEN 530* 531* IJOB=2: Choose all intervals containing eigenvalues. 532* 533* Insure that N(w) is monotone 534* 535 ITMP1 = MIN( NAB( JI, 2 ), 536 $ MAX( NAB( JI, 1 ), ITMP1 ) ) 537* 538* Update the Queue -- add intervals if both halves 539* contain eigenvalues. 540* 541 IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN 542* 543* No eigenvalue in the upper interval: 544* just use the lower interval. 545* 546 AB( JI, 2 ) = TMP1 547* 548 ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN 549* 550* No eigenvalue in the lower interval: 551* just use the upper interval. 552* 553 AB( JI, 1 ) = TMP1 554 ELSE IF( KLNEW.LT.MMAX ) THEN 555* 556* Eigenvalue in both intervals -- add upper to queue. 557* 558 KLNEW = KLNEW + 1 559 AB( KLNEW, 2 ) = AB( JI, 2 ) 560 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 561 AB( KLNEW, 1 ) = TMP1 562 NAB( KLNEW, 1 ) = ITMP1 563 AB( JI, 2 ) = TMP1 564 NAB( JI, 2 ) = ITMP1 565 ELSE 566 INFO = MMAX + 1 567 RETURN 568 END IF 569 ELSE 570* 571* IJOB=3: Binary search. Keep only the interval 572* containing w s.t. N(w) = NVAL 573* 574 IF( ITMP1.LE.NVAL( JI ) ) THEN 575 AB( JI, 1 ) = TMP1 576 NAB( JI, 1 ) = ITMP1 577 END IF 578 IF( ITMP1.GE.NVAL( JI ) ) THEN 579 AB( JI, 2 ) = TMP1 580 NAB( JI, 2 ) = ITMP1 581 END IF 582 END IF 583 100 CONTINUE 584 KL = KLNEW 585* 586 END IF 587* 588* Check for convergence 589* 590 KFNEW = KF 591 DO 110 JI = KF, KL 592 TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) ) 593 TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) ) 594 IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR. 595 $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN 596* 597* Converged -- Swap with position KFNEW, 598* then increment KFNEW 599* 600 IF( JI.GT.KFNEW ) THEN 601 TMP1 = AB( JI, 1 ) 602 TMP2 = AB( JI, 2 ) 603 ITMP1 = NAB( JI, 1 ) 604 ITMP2 = NAB( JI, 2 ) 605 AB( JI, 1 ) = AB( KFNEW, 1 ) 606 AB( JI, 2 ) = AB( KFNEW, 2 ) 607 NAB( JI, 1 ) = NAB( KFNEW, 1 ) 608 NAB( JI, 2 ) = NAB( KFNEW, 2 ) 609 AB( KFNEW, 1 ) = TMP1 610 AB( KFNEW, 2 ) = TMP2 611 NAB( KFNEW, 1 ) = ITMP1 612 NAB( KFNEW, 2 ) = ITMP2 613 IF( IJOB.EQ.3 ) THEN 614 ITMP1 = NVAL( JI ) 615 NVAL( JI ) = NVAL( KFNEW ) 616 NVAL( KFNEW ) = ITMP1 617 END IF 618 END IF 619 KFNEW = KFNEW + 1 620 END IF 621 110 CONTINUE 622 KF = KFNEW 623* 624* Choose Midpoints 625* 626 DO 120 JI = KF, KL 627 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 628 120 CONTINUE 629* 630* If no more intervals to refine, quit. 631* 632 IF( KF.GT.KL ) 633 $ GO TO 140 634 130 CONTINUE 635* 636* Converged 637* 638 140 CONTINUE 639 INFO = MAX( KL+1-KF, 0 ) 640 MOUT = KL 641* 642 RETURN 643* 644* End of DLAEBZ 645* 646 END 647