1*> \brief \b SLAEBZ 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 SLAEBZ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaebz.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaebz.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaebz.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLAEBZ( 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* REAL ABSTOL, PIVMIN, RELTOL 28* .. 29* .. Array Arguments .. 30* INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 31* REAL AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 32* $ WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> SLAEBZ 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 SLAEBZ 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 SLAEBZ 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (N) 167*> The diagonal elements of the tridiagonal matrix T. 168*> \endverbatim 169*> 170*> \param[in] E 171*> \verbatim 172*> E is REAL 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 REAL 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 REAL 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 REAL 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 SLAEBZ 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 SLAEBZ is called. 246*> \endverbatim 247*> 248*> \param[out] WORK 249*> \verbatim 250*> WORK is REAL 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*> \date September 2012 277* 278*> \ingroup auxOTHERauxiliary 279* 280*> \par Further Details: 281* ===================== 282*> 283*> \verbatim 284*> 285*> This routine is intended to be called only by other LAPACK 286*> routines, thus the interface is less user-friendly. It is intended 287*> for two purposes: 288*> 289*> (a) finding eigenvalues. In this case, SLAEBZ should have one or 290*> more initial intervals set up in AB, and SLAEBZ should be called 291*> with IJOB=1. This sets up NAB, and also counts the eigenvalues. 292*> Intervals with no eigenvalues would usually be thrown out at 293*> this point. Also, if not all the eigenvalues in an interval i 294*> are desired, NAB(i,1) can be increased or NAB(i,2) decreased. 295*> For example, set NAB(i,1)=NAB(i,2)-1 to get the largest 296*> eigenvalue. SLAEBZ is then called with IJOB=2 and MMAX 297*> no smaller than the value of MOUT returned by the call with 298*> IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 299*> through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the 300*> tolerance specified by ABSTOL and RELTOL. 301*> 302*> (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). 303*> In this case, start with a Gershgorin interval (a,b). Set up 304*> AB to contain 2 search intervals, both initially (a,b). One 305*> NVAL element should contain f-1 and the other should contain l 306*> , while C should contain a and b, resp. NAB(i,1) should be -1 307*> and NAB(i,2) should be N+1, to flag an error if the desired 308*> interval does not lie in (a,b). SLAEBZ is then called with 309*> IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- 310*> j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while 311*> if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r 312*> >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and 313*> N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and 314*> w(l-r)=...=w(l+k) are handled similarly. 315*> \endverbatim 316*> 317* ===================================================================== 318 SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, 319 $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, 320 $ NAB, WORK, IWORK, INFO ) 321* 322* -- LAPACK auxiliary routine (version 3.4.2) -- 323* -- LAPACK is a software package provided by Univ. of Tennessee, -- 324* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 325* September 2012 326* 327* .. Scalar Arguments .. 328 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX 329 REAL ABSTOL, PIVMIN, RELTOL 330* .. 331* .. Array Arguments .. 332 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 333 REAL AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 334 $ WORK( * ) 335* .. 336* 337* ===================================================================== 338* 339* .. Parameters .. 340 REAL ZERO, TWO, HALF 341 PARAMETER ( ZERO = 0.0E0, TWO = 2.0E0, 342 $ HALF = 1.0E0 / TWO ) 343* .. 344* .. Local Scalars .. 345 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, 346 $ KLNEW 347 REAL TMP1, TMP2 348* .. 349* .. Intrinsic Functions .. 350 INTRINSIC ABS, MAX, MIN 351* .. 352* .. Executable Statements .. 353* 354* Check for Errors 355* 356 INFO = 0 357 IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN 358 INFO = -1 359 RETURN 360 END IF 361* 362* Initialize NAB 363* 364 IF( IJOB.EQ.1 ) THEN 365* 366* Compute the number of eigenvalues in the initial intervals. 367* 368 MOUT = 0 369 DO 30 JI = 1, MINP 370 DO 20 JP = 1, 2 371 TMP1 = D( 1 ) - AB( JI, JP ) 372 IF( ABS( TMP1 ).LT.PIVMIN ) 373 $ TMP1 = -PIVMIN 374 NAB( JI, JP ) = 0 375 IF( TMP1.LE.ZERO ) 376 $ NAB( JI, JP ) = 1 377* 378 DO 10 J = 2, N 379 TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) 380 IF( ABS( TMP1 ).LT.PIVMIN ) 381 $ TMP1 = -PIVMIN 382 IF( TMP1.LE.ZERO ) 383 $ NAB( JI, JP ) = NAB( JI, JP ) + 1 384 10 CONTINUE 385 20 CONTINUE 386 MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) 387 30 CONTINUE 388 RETURN 389 END IF 390* 391* Initialize for loop 392* 393* KF and KL have the following meaning: 394* Intervals 1,...,KF-1 have converged. 395* Intervals KF,...,KL still need to be refined. 396* 397 KF = 1 398 KL = MINP 399* 400* If IJOB=2, initialize C. 401* If IJOB=3, use the user-supplied starting point. 402* 403 IF( IJOB.EQ.2 ) THEN 404 DO 40 JI = 1, MINP 405 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 406 40 CONTINUE 407 END IF 408* 409* Iteration loop 410* 411 DO 130 JIT = 1, NITMAX 412* 413* Loop over intervals 414* 415 IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN 416* 417* Begin of Parallel Version of the loop 418* 419 DO 60 JI = KF, KL 420* 421* Compute N(c), the number of eigenvalues less than c 422* 423 WORK( JI ) = D( 1 ) - C( JI ) 424 IWORK( JI ) = 0 425 IF( WORK( JI ).LE.PIVMIN ) THEN 426 IWORK( JI ) = 1 427 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 428 END IF 429* 430 DO 50 J = 2, N 431 WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) 432 IF( WORK( JI ).LE.PIVMIN ) THEN 433 IWORK( JI ) = IWORK( JI ) + 1 434 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 435 END IF 436 50 CONTINUE 437 60 CONTINUE 438* 439 IF( IJOB.LE.2 ) THEN 440* 441* IJOB=2: Choose all intervals containing eigenvalues. 442* 443 KLNEW = KL 444 DO 70 JI = KF, KL 445* 446* Insure that N(w) is monotone 447* 448 IWORK( JI ) = MIN( NAB( JI, 2 ), 449 $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) 450* 451* Update the Queue -- add intervals if both halves 452* contain eigenvalues. 453* 454 IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN 455* 456* No eigenvalue in the upper interval: 457* just use the lower interval. 458* 459 AB( JI, 2 ) = C( JI ) 460* 461 ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN 462* 463* No eigenvalue in the lower interval: 464* just use the upper interval. 465* 466 AB( JI, 1 ) = C( JI ) 467 ELSE 468 KLNEW = KLNEW + 1 469 IF( KLNEW.LE.MMAX ) THEN 470* 471* Eigenvalue in both intervals -- add upper to 472* queue. 473* 474 AB( KLNEW, 2 ) = AB( JI, 2 ) 475 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 476 AB( KLNEW, 1 ) = C( JI ) 477 NAB( KLNEW, 1 ) = IWORK( JI ) 478 AB( JI, 2 ) = C( JI ) 479 NAB( JI, 2 ) = IWORK( JI ) 480 ELSE 481 INFO = MMAX + 1 482 END IF 483 END IF 484 70 CONTINUE 485 IF( INFO.NE.0 ) 486 $ RETURN 487 KL = KLNEW 488 ELSE 489* 490* IJOB=3: Binary search. Keep only the interval containing 491* w s.t. N(w) = NVAL 492* 493 DO 80 JI = KF, KL 494 IF( IWORK( JI ).LE.NVAL( JI ) ) THEN 495 AB( JI, 1 ) = C( JI ) 496 NAB( JI, 1 ) = IWORK( JI ) 497 END IF 498 IF( IWORK( JI ).GE.NVAL( JI ) ) THEN 499 AB( JI, 2 ) = C( JI ) 500 NAB( JI, 2 ) = IWORK( JI ) 501 END IF 502 80 CONTINUE 503 END IF 504* 505 ELSE 506* 507* End of Parallel Version of the loop 508* 509* Begin of Serial Version of the loop 510* 511 KLNEW = KL 512 DO 100 JI = KF, KL 513* 514* Compute N(w), the number of eigenvalues less than w 515* 516 TMP1 = C( JI ) 517 TMP2 = D( 1 ) - TMP1 518 ITMP1 = 0 519 IF( TMP2.LE.PIVMIN ) THEN 520 ITMP1 = 1 521 TMP2 = MIN( TMP2, -PIVMIN ) 522 END IF 523* 524 DO 90 J = 2, N 525 TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 526 IF( TMP2.LE.PIVMIN ) THEN 527 ITMP1 = ITMP1 + 1 528 TMP2 = MIN( TMP2, -PIVMIN ) 529 END IF 530 90 CONTINUE 531* 532 IF( IJOB.LE.2 ) THEN 533* 534* IJOB=2: Choose all intervals containing eigenvalues. 535* 536* Insure that N(w) is monotone 537* 538 ITMP1 = MIN( NAB( JI, 2 ), 539 $ MAX( NAB( JI, 1 ), ITMP1 ) ) 540* 541* Update the Queue -- add intervals if both halves 542* contain eigenvalues. 543* 544 IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN 545* 546* No eigenvalue in the upper interval: 547* just use the lower interval. 548* 549 AB( JI, 2 ) = TMP1 550* 551 ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN 552* 553* No eigenvalue in the lower interval: 554* just use the upper interval. 555* 556 AB( JI, 1 ) = TMP1 557 ELSE IF( KLNEW.LT.MMAX ) THEN 558* 559* Eigenvalue in both intervals -- add upper to queue. 560* 561 KLNEW = KLNEW + 1 562 AB( KLNEW, 2 ) = AB( JI, 2 ) 563 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 564 AB( KLNEW, 1 ) = TMP1 565 NAB( KLNEW, 1 ) = ITMP1 566 AB( JI, 2 ) = TMP1 567 NAB( JI, 2 ) = ITMP1 568 ELSE 569 INFO = MMAX + 1 570 RETURN 571 END IF 572 ELSE 573* 574* IJOB=3: Binary search. Keep only the interval 575* containing w s.t. N(w) = NVAL 576* 577 IF( ITMP1.LE.NVAL( JI ) ) THEN 578 AB( JI, 1 ) = TMP1 579 NAB( JI, 1 ) = ITMP1 580 END IF 581 IF( ITMP1.GE.NVAL( JI ) ) THEN 582 AB( JI, 2 ) = TMP1 583 NAB( JI, 2 ) = ITMP1 584 END IF 585 END IF 586 100 CONTINUE 587 KL = KLNEW 588* 589 END IF 590* 591* Check for convergence 592* 593 KFNEW = KF 594 DO 110 JI = KF, KL 595 TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) ) 596 TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) ) 597 IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR. 598 $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN 599* 600* Converged -- Swap with position KFNEW, 601* then increment KFNEW 602* 603 IF( JI.GT.KFNEW ) THEN 604 TMP1 = AB( JI, 1 ) 605 TMP2 = AB( JI, 2 ) 606 ITMP1 = NAB( JI, 1 ) 607 ITMP2 = NAB( JI, 2 ) 608 AB( JI, 1 ) = AB( KFNEW, 1 ) 609 AB( JI, 2 ) = AB( KFNEW, 2 ) 610 NAB( JI, 1 ) = NAB( KFNEW, 1 ) 611 NAB( JI, 2 ) = NAB( KFNEW, 2 ) 612 AB( KFNEW, 1 ) = TMP1 613 AB( KFNEW, 2 ) = TMP2 614 NAB( KFNEW, 1 ) = ITMP1 615 NAB( KFNEW, 2 ) = ITMP2 616 IF( IJOB.EQ.3 ) THEN 617 ITMP1 = NVAL( JI ) 618 NVAL( JI ) = NVAL( KFNEW ) 619 NVAL( KFNEW ) = ITMP1 620 END IF 621 END IF 622 KFNEW = KFNEW + 1 623 END IF 624 110 CONTINUE 625 KF = KFNEW 626* 627* Choose Midpoints 628* 629 DO 120 JI = KF, KL 630 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 631 120 CONTINUE 632* 633* If no more intervals to refine, quit. 634* 635 IF( KF.GT.KL ) 636 $ GO TO 140 637 130 CONTINUE 638* 639* Converged 640* 641 140 CONTINUE 642 INFO = MAX( KL+1-KF, 0 ) 643 MOUT = KL 644* 645 RETURN 646* 647* End of SLAEBZ 648* 649 END 650c $Id$ 651