1*> \brief \b SSTEBZ 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SSTEBZ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstebz.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstebz.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstebz.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, 22* M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, 23* INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER ORDER, RANGE 27* INTEGER IL, INFO, IU, M, N, NSPLIT 28* REAL ABSTOL, VL, VU 29* .. 30* .. Array Arguments .. 31* INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 32* REAL D( * ), E( * ), W( * ), WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> SSTEBZ computes the eigenvalues of a symmetric tridiagonal 42*> matrix T. The user may ask for all eigenvalues, all eigenvalues 43*> in the half-open interval (VL, VU], or the IL-th through IU-th 44*> eigenvalues. 45*> 46*> To avoid overflow, the matrix must be scaled so that its 47*> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest 48*> accuracy, it should not be much smaller than that. 49*> 50*> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 51*> Matrix", Report CS41, Computer Science Dept., Stanford 52*> University, July 21, 1966. 53*> \endverbatim 54* 55* Arguments: 56* ========== 57* 58*> \param[in] RANGE 59*> \verbatim 60*> RANGE is CHARACTER*1 61*> = 'A': ("All") all eigenvalues will be found. 62*> = 'V': ("Value") all eigenvalues in the half-open interval 63*> (VL, VU] will be found. 64*> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the 65*> entire matrix) will be found. 66*> \endverbatim 67*> 68*> \param[in] ORDER 69*> \verbatim 70*> ORDER is CHARACTER*1 71*> = 'B': ("By Block") the eigenvalues will be grouped by 72*> split-off block (see IBLOCK, ISPLIT) and 73*> ordered from smallest to largest within 74*> the block. 75*> = 'E': ("Entire matrix") 76*> the eigenvalues for the entire matrix 77*> will be ordered from smallest to 78*> largest. 79*> \endverbatim 80*> 81*> \param[in] N 82*> \verbatim 83*> N is INTEGER 84*> The order of the tridiagonal matrix T. N >= 0. 85*> \endverbatim 86*> 87*> \param[in] VL 88*> \verbatim 89*> VL is REAL 90*> 91*> If RANGE='V', the lower bound of the interval to 92*> be searched for eigenvalues. Eigenvalues less than or equal 93*> to VL, or greater than VU, will not be returned. VL < VU. 94*> Not referenced if RANGE = 'A' or 'I'. 95*> \endverbatim 96*> 97*> \param[in] VU 98*> \verbatim 99*> VU is REAL 100*> 101*> If RANGE='V', the upper bound of the interval to 102*> be searched for eigenvalues. Eigenvalues less than or equal 103*> to VL, or greater than VU, will not be returned. VL < VU. 104*> Not referenced if RANGE = 'A' or 'I'. 105*> \endverbatim 106*> 107*> \param[in] IL 108*> \verbatim 109*> IL is INTEGER 110*> 111*> If RANGE='I', the index of the 112*> smallest eigenvalue to be returned. 113*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 114*> Not referenced if RANGE = 'A' or 'V'. 115*> \endverbatim 116*> 117*> \param[in] IU 118*> \verbatim 119*> IU is INTEGER 120*> 121*> If RANGE='I', the index of the 122*> largest eigenvalue to be returned. 123*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 124*> Not referenced if RANGE = 'A' or 'V'. 125*> \endverbatim 126*> 127*> \param[in] ABSTOL 128*> \verbatim 129*> ABSTOL is REAL 130*> The absolute tolerance for the eigenvalues. An eigenvalue 131*> (or cluster) is considered to be located if it has been 132*> determined to lie in an interval whose width is ABSTOL or 133*> less. If ABSTOL is less than or equal to zero, then ULP*|T| 134*> will be used, where |T| means the 1-norm of T. 135*> 136*> Eigenvalues will be computed most accurately when ABSTOL is 137*> set to twice the underflow threshold 2*SLAMCH('S'), not zero. 138*> \endverbatim 139*> 140*> \param[in] D 141*> \verbatim 142*> D is REAL array, dimension (N) 143*> The n diagonal elements of the tridiagonal matrix T. 144*> \endverbatim 145*> 146*> \param[in] E 147*> \verbatim 148*> E is REAL array, dimension (N-1) 149*> The (n-1) off-diagonal elements of the tridiagonal matrix T. 150*> \endverbatim 151*> 152*> \param[out] M 153*> \verbatim 154*> M is INTEGER 155*> The actual number of eigenvalues found. 0 <= M <= N. 156*> (See also the description of INFO=2,3.) 157*> \endverbatim 158*> 159*> \param[out] NSPLIT 160*> \verbatim 161*> NSPLIT is INTEGER 162*> The number of diagonal blocks in the matrix T. 163*> 1 <= NSPLIT <= N. 164*> \endverbatim 165*> 166*> \param[out] W 167*> \verbatim 168*> W is REAL array, dimension (N) 169*> On exit, the first M elements of W will contain the 170*> eigenvalues. (SSTEBZ may use the remaining N-M elements as 171*> workspace.) 172*> \endverbatim 173*> 174*> \param[out] IBLOCK 175*> \verbatim 176*> IBLOCK is INTEGER array, dimension (N) 177*> At each row/column j where E(j) is zero or small, the 178*> matrix T is considered to split into a block diagonal 179*> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which 180*> block (from 1 to the number of blocks) the eigenvalue W(i) 181*> belongs. (SSTEBZ may use the remaining N-M elements as 182*> workspace.) 183*> \endverbatim 184*> 185*> \param[out] ISPLIT 186*> \verbatim 187*> ISPLIT is INTEGER array, dimension (N) 188*> The splitting points, at which T breaks up into submatrices. 189*> The first submatrix consists of rows/columns 1 to ISPLIT(1), 190*> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), 191*> etc., and the NSPLIT-th consists of rows/columns 192*> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 193*> (Only the first NSPLIT elements will actually be used, but 194*> since the user cannot know a priori what value NSPLIT will 195*> have, N words must be reserved for ISPLIT.) 196*> \endverbatim 197*> 198*> \param[out] WORK 199*> \verbatim 200*> WORK is REAL array, dimension (4*N) 201*> \endverbatim 202*> 203*> \param[out] IWORK 204*> \verbatim 205*> IWORK is INTEGER array, dimension (3*N) 206*> \endverbatim 207*> 208*> \param[out] INFO 209*> \verbatim 210*> INFO is INTEGER 211*> = 0: successful exit 212*> < 0: if INFO = -i, the i-th argument had an illegal value 213*> > 0: some or all of the eigenvalues failed to converge or 214*> were not computed: 215*> =1 or 3: Bisection failed to converge for some 216*> eigenvalues; these eigenvalues are flagged by a 217*> negative block number. The effect is that the 218*> eigenvalues may not be as accurate as the 219*> absolute and relative tolerances. This is 220*> generally caused by unexpectedly inaccurate 221*> arithmetic. 222*> =2 or 3: RANGE='I' only: Not all of the eigenvalues 223*> IL:IU were found. 224*> Effect: M < IU+1-IL 225*> Cause: non-monotonic arithmetic, causing the 226*> Sturm sequence to be non-monotonic. 227*> Cure: recalculate, using RANGE='A', and pick 228*> out eigenvalues IL:IU. In some cases, 229*> increasing the PARAMETER "FUDGE" may 230*> make things work. 231*> = 4: RANGE='I', and the Gershgorin interval 232*> initially used was too small. No eigenvalues 233*> were computed. 234*> Probable cause: your machine has sloppy 235*> floating-point arithmetic. 236*> Cure: Increase the PARAMETER "FUDGE", 237*> recompile, and try again. 238*> \endverbatim 239* 240*> \par Internal Parameters: 241* ========================= 242*> 243*> \verbatim 244*> RELFAC REAL, default = 2.0e0 245*> The relative tolerance. An interval (a,b] lies within 246*> "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), 247*> where "ulp" is the machine precision (distance from 1 to 248*> the next larger floating point number.) 249*> 250*> FUDGE REAL, default = 2 251*> A "fudge factor" to widen the Gershgorin intervals. Ideally, 252*> a value of 1 should work, but on machines with sloppy 253*> arithmetic, this needs to be larger. The default for 254*> publicly released versions should be large enough to handle 255*> the worst machine around. Note that this has no effect 256*> on accuracy of the solution. 257*> \endverbatim 258* 259* Authors: 260* ======== 261* 262*> \author Univ. of Tennessee 263*> \author Univ. of California Berkeley 264*> \author Univ. of Colorado Denver 265*> \author NAG Ltd. 266* 267*> \ingroup auxOTHERcomputational 268* 269* ===================================================================== 270 SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, 271 $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, 272 $ INFO ) 273* 274* -- LAPACK computational routine -- 275* -- LAPACK is a software package provided by Univ. of Tennessee, -- 276* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 277* 278* .. Scalar Arguments .. 279 CHARACTER ORDER, RANGE 280 INTEGER IL, INFO, IU, M, N, NSPLIT 281 REAL ABSTOL, VL, VU 282* .. 283* .. Array Arguments .. 284 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 285 REAL D( * ), E( * ), W( * ), WORK( * ) 286* .. 287* 288* ===================================================================== 289* 290* .. Parameters .. 291 REAL ZERO, ONE, TWO, HALF 292 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 293 $ HALF = 1.0E0 / TWO ) 294 REAL FUDGE, RELFAC 295 PARAMETER ( FUDGE = 2.1E0, RELFAC = 2.0E0 ) 296* .. 297* .. Local Scalars .. 298 LOGICAL NCNVRG, TOOFEW 299 INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO, 300 $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX, 301 $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL, 302 $ NWU 303 REAL ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN, 304 $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL 305* .. 306* .. Local Arrays .. 307 INTEGER IDUMMA( 1 ) 308* .. 309* .. External Functions .. 310 LOGICAL LSAME 311 INTEGER ILAENV 312 REAL SLAMCH 313 EXTERNAL LSAME, ILAENV, SLAMCH 314* .. 315* .. External Subroutines .. 316 EXTERNAL SLAEBZ, XERBLA 317* .. 318* .. Intrinsic Functions .. 319 INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT 320* .. 321* .. Executable Statements .. 322* 323 INFO = 0 324* 325* Decode RANGE 326* 327 IF( LSAME( RANGE, 'A' ) ) THEN 328 IRANGE = 1 329 ELSE IF( LSAME( RANGE, 'V' ) ) THEN 330 IRANGE = 2 331 ELSE IF( LSAME( RANGE, 'I' ) ) THEN 332 IRANGE = 3 333 ELSE 334 IRANGE = 0 335 END IF 336* 337* Decode ORDER 338* 339 IF( LSAME( ORDER, 'B' ) ) THEN 340 IORDER = 2 341 ELSE IF( LSAME( ORDER, 'E' ) ) THEN 342 IORDER = 1 343 ELSE 344 IORDER = 0 345 END IF 346* 347* Check for Errors 348* 349 IF( IRANGE.LE.0 ) THEN 350 INFO = -1 351 ELSE IF( IORDER.LE.0 ) THEN 352 INFO = -2 353 ELSE IF( N.LT.0 ) THEN 354 INFO = -3 355 ELSE IF( IRANGE.EQ.2 ) THEN 356 IF( VL.GE.VU ) INFO = -5 357 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) 358 $ THEN 359 INFO = -6 360 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) 361 $ THEN 362 INFO = -7 363 END IF 364* 365 IF( INFO.NE.0 ) THEN 366 CALL XERBLA( 'SSTEBZ', -INFO ) 367 RETURN 368 END IF 369* 370* Initialize error flags 371* 372 INFO = 0 373 NCNVRG = .FALSE. 374 TOOFEW = .FALSE. 375* 376* Quick return if possible 377* 378 M = 0 379 IF( N.EQ.0 ) 380 $ RETURN 381* 382* Simplifications: 383* 384 IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N ) 385 $ IRANGE = 1 386* 387* Get machine constants 388* NB is the minimum vector length for vector bisection, or 0 389* if only scalar is to be done. 390* 391 SAFEMN = SLAMCH( 'S' ) 392 ULP = SLAMCH( 'P' ) 393 RTOLI = ULP*RELFAC 394 NB = ILAENV( 1, 'SSTEBZ', ' ', N, -1, -1, -1 ) 395 IF( NB.LE.1 ) 396 $ NB = 0 397* 398* Special Case when N=1 399* 400 IF( N.EQ.1 ) THEN 401 NSPLIT = 1 402 ISPLIT( 1 ) = 1 403 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN 404 M = 0 405 ELSE 406 W( 1 ) = D( 1 ) 407 IBLOCK( 1 ) = 1 408 M = 1 409 END IF 410 RETURN 411 END IF 412* 413* Compute Splitting Points 414* 415 NSPLIT = 1 416 WORK( N ) = ZERO 417 PIVMIN = ONE 418* 419 DO 10 J = 2, N 420 TMP1 = E( J-1 )**2 421 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN 422 ISPLIT( NSPLIT ) = J - 1 423 NSPLIT = NSPLIT + 1 424 WORK( J-1 ) = ZERO 425 ELSE 426 WORK( J-1 ) = TMP1 427 PIVMIN = MAX( PIVMIN, TMP1 ) 428 END IF 429 10 CONTINUE 430 ISPLIT( NSPLIT ) = N 431 PIVMIN = PIVMIN*SAFEMN 432* 433* Compute Interval and ATOLI 434* 435 IF( IRANGE.EQ.3 ) THEN 436* 437* RANGE='I': Compute the interval containing eigenvalues 438* IL through IU. 439* 440* Compute Gershgorin interval for entire (split) matrix 441* and use it as the initial interval 442* 443 GU = D( 1 ) 444 GL = D( 1 ) 445 TMP1 = ZERO 446* 447 DO 20 J = 1, N - 1 448 TMP2 = SQRT( WORK( J ) ) 449 GU = MAX( GU, D( J )+TMP1+TMP2 ) 450 GL = MIN( GL, D( J )-TMP1-TMP2 ) 451 TMP1 = TMP2 452 20 CONTINUE 453* 454 GU = MAX( GU, D( N )+TMP1 ) 455 GL = MIN( GL, D( N )-TMP1 ) 456 TNORM = MAX( ABS( GL ), ABS( GU ) ) 457 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN 458 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN 459* 460* Compute Iteration parameters 461* 462 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 463 $ LOG( TWO ) ) + 2 464 IF( ABSTOL.LE.ZERO ) THEN 465 ATOLI = ULP*TNORM 466 ELSE 467 ATOLI = ABSTOL 468 END IF 469* 470 WORK( N+1 ) = GL 471 WORK( N+2 ) = GL 472 WORK( N+3 ) = GU 473 WORK( N+4 ) = GU 474 WORK( N+5 ) = GL 475 WORK( N+6 ) = GU 476 IWORK( 1 ) = -1 477 IWORK( 2 ) = -1 478 IWORK( 3 ) = N + 1 479 IWORK( 4 ) = N + 1 480 IWORK( 5 ) = IL - 1 481 IWORK( 6 ) = IU 482* 483 CALL SLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, 484 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 485 $ IWORK, W, IBLOCK, IINFO ) 486* 487 IF( IWORK( 6 ).EQ.IU ) THEN 488 WL = WORK( N+1 ) 489 WLU = WORK( N+3 ) 490 NWL = IWORK( 1 ) 491 WU = WORK( N+4 ) 492 WUL = WORK( N+2 ) 493 NWU = IWORK( 4 ) 494 ELSE 495 WL = WORK( N+2 ) 496 WLU = WORK( N+4 ) 497 NWL = IWORK( 2 ) 498 WU = WORK( N+3 ) 499 WUL = WORK( N+1 ) 500 NWU = IWORK( 3 ) 501 END IF 502* 503 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 504 INFO = 4 505 RETURN 506 END IF 507 ELSE 508* 509* RANGE='A' or 'V' -- Set ATOLI 510* 511 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), 512 $ ABS( D( N ) )+ABS( E( N-1 ) ) ) 513* 514 DO 30 J = 2, N - 1 515 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ 516 $ ABS( E( J ) ) ) 517 30 CONTINUE 518* 519 IF( ABSTOL.LE.ZERO ) THEN 520 ATOLI = ULP*TNORM 521 ELSE 522 ATOLI = ABSTOL 523 END IF 524* 525 IF( IRANGE.EQ.2 ) THEN 526 WL = VL 527 WU = VU 528 ELSE 529 WL = ZERO 530 WU = ZERO 531 END IF 532 END IF 533* 534* Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 535* NWL accumulates the number of eigenvalues .le. WL, 536* NWU accumulates the number of eigenvalues .le. WU 537* 538 M = 0 539 IEND = 0 540 INFO = 0 541 NWL = 0 542 NWU = 0 543* 544 DO 70 JB = 1, NSPLIT 545 IOFF = IEND 546 IBEGIN = IOFF + 1 547 IEND = ISPLIT( JB ) 548 IN = IEND - IOFF 549* 550 IF( IN.EQ.1 ) THEN 551* 552* Special Case -- IN=1 553* 554 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) 555 $ NWL = NWL + 1 556 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) 557 $ NWU = NWU + 1 558 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. 559 $ D( IBEGIN )-PIVMIN ) ) THEN 560 M = M + 1 561 W( M ) = D( IBEGIN ) 562 IBLOCK( M ) = JB 563 END IF 564 ELSE 565* 566* General Case -- IN > 1 567* 568* Compute Gershgorin Interval 569* and use it as the initial interval 570* 571 GU = D( IBEGIN ) 572 GL = D( IBEGIN ) 573 TMP1 = ZERO 574* 575 DO 40 J = IBEGIN, IEND - 1 576 TMP2 = ABS( E( J ) ) 577 GU = MAX( GU, D( J )+TMP1+TMP2 ) 578 GL = MIN( GL, D( J )-TMP1-TMP2 ) 579 TMP1 = TMP2 580 40 CONTINUE 581* 582 GU = MAX( GU, D( IEND )+TMP1 ) 583 GL = MIN( GL, D( IEND )-TMP1 ) 584 BNORM = MAX( ABS( GL ), ABS( GU ) ) 585 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN 586 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN 587* 588* Compute ATOLI for the current submatrix 589* 590 IF( ABSTOL.LE.ZERO ) THEN 591 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) ) 592 ELSE 593 ATOLI = ABSTOL 594 END IF 595* 596 IF( IRANGE.GT.1 ) THEN 597 IF( GU.LT.WL ) THEN 598 NWL = NWL + IN 599 NWU = NWU + IN 600 GO TO 70 601 END IF 602 GL = MAX( GL, WL ) 603 GU = MIN( GU, WU ) 604 IF( GL.GE.GU ) 605 $ GO TO 70 606 END IF 607* 608* Set Up Initial Interval 609* 610 WORK( N+1 ) = GL 611 WORK( N+IN+1 ) = GU 612 CALL SLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 613 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 614 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 615 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 616* 617 NWL = NWL + IWORK( 1 ) 618 NWU = NWU + IWORK( IN+1 ) 619 IWOFF = M - IWORK( 1 ) 620* 621* Compute Eigenvalues 622* 623 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 624 $ LOG( TWO ) ) + 2 625 CALL SLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 626 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 627 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 628 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 629* 630* Copy Eigenvalues Into W and IBLOCK 631* Use -JB for block number for unconverged eigenvalues. 632* 633 DO 60 J = 1, IOUT 634 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 635* 636* Flag non-convergence. 637* 638 IF( J.GT.IOUT-IINFO ) THEN 639 NCNVRG = .TRUE. 640 IB = -JB 641 ELSE 642 IB = JB 643 END IF 644 DO 50 JE = IWORK( J ) + 1 + IWOFF, 645 $ IWORK( J+IN ) + IWOFF 646 W( JE ) = TMP1 647 IBLOCK( JE ) = IB 648 50 CONTINUE 649 60 CONTINUE 650* 651 M = M + IM 652 END IF 653 70 CONTINUE 654* 655* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 656* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 657* 658 IF( IRANGE.EQ.3 ) THEN 659 IM = 0 660 IDISCL = IL - 1 - NWL 661 IDISCU = NWU - IU 662* 663 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 664 DO 80 JE = 1, M 665 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 666 IDISCL = IDISCL - 1 667 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 668 IDISCU = IDISCU - 1 669 ELSE 670 IM = IM + 1 671 W( IM ) = W( JE ) 672 IBLOCK( IM ) = IBLOCK( JE ) 673 END IF 674 80 CONTINUE 675 M = IM 676 END IF 677 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 678* 679* Code to deal with effects of bad arithmetic: 680* Some low eigenvalues to be discarded are not in (WL,WLU], 681* or high eigenvalues to be discarded are not in (WUL,WU] 682* so just kill off the smallest IDISCL/largest IDISCU 683* eigenvalues, by simply finding the smallest/largest 684* eigenvalue(s). 685* 686* (If N(w) is monotone non-decreasing, this should never 687* happen.) 688* 689 IF( IDISCL.GT.0 ) THEN 690 WKILL = WU 691 DO 100 JDISC = 1, IDISCL 692 IW = 0 693 DO 90 JE = 1, M 694 IF( IBLOCK( JE ).NE.0 .AND. 695 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 696 IW = JE 697 WKILL = W( JE ) 698 END IF 699 90 CONTINUE 700 IBLOCK( IW ) = 0 701 100 CONTINUE 702 END IF 703 IF( IDISCU.GT.0 ) THEN 704* 705 WKILL = WL 706 DO 120 JDISC = 1, IDISCU 707 IW = 0 708 DO 110 JE = 1, M 709 IF( IBLOCK( JE ).NE.0 .AND. 710 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN 711 IW = JE 712 WKILL = W( JE ) 713 END IF 714 110 CONTINUE 715 IBLOCK( IW ) = 0 716 120 CONTINUE 717 END IF 718 IM = 0 719 DO 130 JE = 1, M 720 IF( IBLOCK( JE ).NE.0 ) THEN 721 IM = IM + 1 722 W( IM ) = W( JE ) 723 IBLOCK( IM ) = IBLOCK( JE ) 724 END IF 725 130 CONTINUE 726 M = IM 727 END IF 728 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 729 TOOFEW = .TRUE. 730 END IF 731 END IF 732* 733* If ORDER='B', do nothing -- the eigenvalues are already sorted 734* by block. 735* If ORDER='E', sort the eigenvalues from smallest to largest 736* 737 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN 738 DO 150 JE = 1, M - 1 739 IE = 0 740 TMP1 = W( JE ) 741 DO 140 J = JE + 1, M 742 IF( W( J ).LT.TMP1 ) THEN 743 IE = J 744 TMP1 = W( J ) 745 END IF 746 140 CONTINUE 747* 748 IF( IE.NE.0 ) THEN 749 ITMP1 = IBLOCK( IE ) 750 W( IE ) = W( JE ) 751 IBLOCK( IE ) = IBLOCK( JE ) 752 W( JE ) = TMP1 753 IBLOCK( JE ) = ITMP1 754 END IF 755 150 CONTINUE 756 END IF 757* 758 INFO = 0 759 IF( NCNVRG ) 760 $ INFO = INFO + 1 761 IF( TOOFEW ) 762 $ INFO = INFO + 2 763 RETURN 764* 765* End of SSTEBZ 766* 767 END 768