1*> \brief \b DSTEBZ 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DSTEBZ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstebz.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstebz.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstebz.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSTEBZ( 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* DOUBLE PRECISION ABSTOL, VL, VU 29* .. 30* .. Array Arguments .. 31* INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 32* DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> DSTEBZ 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*DLAMCH('S'), not zero. 138*> \endverbatim 139*> 140*> \param[in] D 141*> \verbatim 142*> D is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) 169*> On exit, the first M elements of W will contain the 170*> eigenvalues. (DSTEBZ 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. (DSTEBZ 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 DOUBLE PRECISION 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 DOUBLE PRECISION, 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 DOUBLE PRECISION, 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 DSTEBZ( 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 DOUBLE PRECISION ABSTOL, VL, VU 282* .. 283* .. Array Arguments .. 284 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 285 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 286* .. 287* 288* ===================================================================== 289* 290* .. Parameters .. 291 DOUBLE PRECISION ZERO, ONE, TWO, HALF 292 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 293 $ HALF = 1.0D0 / TWO ) 294 DOUBLE PRECISION FUDGE, RELFAC 295 PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 ) 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 313 EXTERNAL LSAME, ILAENV, DLAMCH 314* .. 315* .. External Subroutines .. 316 EXTERNAL DLAEBZ, 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 ) 357 $ INFO = -5 358 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) 359 $ THEN 360 INFO = -6 361 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) 362 $ THEN 363 INFO = -7 364 END IF 365* 366 IF( INFO.NE.0 ) THEN 367 CALL XERBLA( 'DSTEBZ', -INFO ) 368 RETURN 369 END IF 370* 371* Initialize error flags 372* 373 INFO = 0 374 NCNVRG = .FALSE. 375 TOOFEW = .FALSE. 376* 377* Quick return if possible 378* 379 M = 0 380 IF( N.EQ.0 ) 381 $ RETURN 382* 383* Simplifications: 384* 385 IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N ) 386 $ IRANGE = 1 387* 388* Get machine constants 389* NB is the minimum vector length for vector bisection, or 0 390* if only scalar is to be done. 391* 392 SAFEMN = DLAMCH( 'S' ) 393 ULP = DLAMCH( 'P' ) 394 RTOLI = ULP*RELFAC 395 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 ) 396 IF( NB.LE.1 ) 397 $ NB = 0 398* 399* Special Case when N=1 400* 401 IF( N.EQ.1 ) THEN 402 NSPLIT = 1 403 ISPLIT( 1 ) = 1 404 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN 405 M = 0 406 ELSE 407 W( 1 ) = D( 1 ) 408 IBLOCK( 1 ) = 1 409 M = 1 410 END IF 411 RETURN 412 END IF 413* 414* Compute Splitting Points 415* 416 NSPLIT = 1 417 WORK( N ) = ZERO 418 PIVMIN = ONE 419* 420 DO 10 J = 2, N 421 TMP1 = E( J-1 )**2 422 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN 423 ISPLIT( NSPLIT ) = J - 1 424 NSPLIT = NSPLIT + 1 425 WORK( J-1 ) = ZERO 426 ELSE 427 WORK( J-1 ) = TMP1 428 PIVMIN = MAX( PIVMIN, TMP1 ) 429 END IF 430 10 CONTINUE 431 ISPLIT( NSPLIT ) = N 432 PIVMIN = PIVMIN*SAFEMN 433* 434* Compute Interval and ATOLI 435* 436 IF( IRANGE.EQ.3 ) THEN 437* 438* RANGE='I': Compute the interval containing eigenvalues 439* IL through IU. 440* 441* Compute Gershgorin interval for entire (split) matrix 442* and use it as the initial interval 443* 444 GU = D( 1 ) 445 GL = D( 1 ) 446 TMP1 = ZERO 447* 448 DO 20 J = 1, N - 1 449 TMP2 = SQRT( WORK( J ) ) 450 GU = MAX( GU, D( J )+TMP1+TMP2 ) 451 GL = MIN( GL, D( J )-TMP1-TMP2 ) 452 TMP1 = TMP2 453 20 CONTINUE 454* 455 GU = MAX( GU, D( N )+TMP1 ) 456 GL = MIN( GL, D( N )-TMP1 ) 457 TNORM = MAX( ABS( GL ), ABS( GU ) ) 458 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN 459 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN 460* 461* Compute Iteration parameters 462* 463 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 464 $ LOG( TWO ) ) + 2 465 IF( ABSTOL.LE.ZERO ) THEN 466 ATOLI = ULP*TNORM 467 ELSE 468 ATOLI = ABSTOL 469 END IF 470* 471 WORK( N+1 ) = GL 472 WORK( N+2 ) = GL 473 WORK( N+3 ) = GU 474 WORK( N+4 ) = GU 475 WORK( N+5 ) = GL 476 WORK( N+6 ) = GU 477 IWORK( 1 ) = -1 478 IWORK( 2 ) = -1 479 IWORK( 3 ) = N + 1 480 IWORK( 4 ) = N + 1 481 IWORK( 5 ) = IL - 1 482 IWORK( 6 ) = IU 483* 484 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, 485 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 486 $ IWORK, W, IBLOCK, IINFO ) 487* 488 IF( IWORK( 6 ).EQ.IU ) THEN 489 WL = WORK( N+1 ) 490 WLU = WORK( N+3 ) 491 NWL = IWORK( 1 ) 492 WU = WORK( N+4 ) 493 WUL = WORK( N+2 ) 494 NWU = IWORK( 4 ) 495 ELSE 496 WL = WORK( N+2 ) 497 WLU = WORK( N+4 ) 498 NWL = IWORK( 2 ) 499 WU = WORK( N+3 ) 500 WUL = WORK( N+1 ) 501 NWU = IWORK( 3 ) 502 END IF 503* 504 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 505 INFO = 4 506 RETURN 507 END IF 508 ELSE 509* 510* RANGE='A' or 'V' -- Set ATOLI 511* 512 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), 513 $ ABS( D( N ) )+ABS( E( N-1 ) ) ) 514* 515 DO 30 J = 2, N - 1 516 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ 517 $ ABS( E( J ) ) ) 518 30 CONTINUE 519* 520 IF( ABSTOL.LE.ZERO ) THEN 521 ATOLI = ULP*TNORM 522 ELSE 523 ATOLI = ABSTOL 524 END IF 525* 526 IF( IRANGE.EQ.2 ) THEN 527 WL = VL 528 WU = VU 529 ELSE 530 WL = ZERO 531 WU = ZERO 532 END IF 533 END IF 534* 535* Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 536* NWL accumulates the number of eigenvalues .le. WL, 537* NWU accumulates the number of eigenvalues .le. WU 538* 539 M = 0 540 IEND = 0 541 INFO = 0 542 NWL = 0 543 NWU = 0 544* 545 DO 70 JB = 1, NSPLIT 546 IOFF = IEND 547 IBEGIN = IOFF + 1 548 IEND = ISPLIT( JB ) 549 IN = IEND - IOFF 550* 551 IF( IN.EQ.1 ) THEN 552* 553* Special Case -- IN=1 554* 555 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) 556 $ NWL = NWL + 1 557 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) 558 $ NWU = NWU + 1 559 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. 560 $ D( IBEGIN )-PIVMIN ) ) THEN 561 M = M + 1 562 W( M ) = D( IBEGIN ) 563 IBLOCK( M ) = JB 564 END IF 565 ELSE 566* 567* General Case -- IN > 1 568* 569* Compute Gershgorin Interval 570* and use it as the initial interval 571* 572 GU = D( IBEGIN ) 573 GL = D( IBEGIN ) 574 TMP1 = ZERO 575* 576 DO 40 J = IBEGIN, IEND - 1 577 TMP2 = ABS( E( J ) ) 578 GU = MAX( GU, D( J )+TMP1+TMP2 ) 579 GL = MIN( GL, D( J )-TMP1-TMP2 ) 580 TMP1 = TMP2 581 40 CONTINUE 582* 583 GU = MAX( GU, D( IEND )+TMP1 ) 584 GL = MIN( GL, D( IEND )-TMP1 ) 585 BNORM = MAX( ABS( GL ), ABS( GU ) ) 586 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN 587 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN 588* 589* Compute ATOLI for the current submatrix 590* 591 IF( ABSTOL.LE.ZERO ) THEN 592 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) ) 593 ELSE 594 ATOLI = ABSTOL 595 END IF 596* 597 IF( IRANGE.GT.1 ) THEN 598 IF( GU.LT.WL ) THEN 599 NWL = NWL + IN 600 NWU = NWU + IN 601 GO TO 70 602 END IF 603 GL = MAX( GL, WL ) 604 GU = MIN( GU, WU ) 605 IF( GL.GE.GU ) 606 $ GO TO 70 607 END IF 608* 609* Set Up Initial Interval 610* 611 WORK( N+1 ) = GL 612 WORK( N+IN+1 ) = GU 613 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 614 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 615 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 616 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 617* 618 NWL = NWL + IWORK( 1 ) 619 NWU = NWU + IWORK( IN+1 ) 620 IWOFF = M - IWORK( 1 ) 621* 622* Compute Eigenvalues 623* 624 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 625 $ LOG( TWO ) ) + 2 626 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 627 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 628 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 629 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 630* 631* Copy Eigenvalues Into W and IBLOCK 632* Use -JB for block number for unconverged eigenvalues. 633* 634 DO 60 J = 1, IOUT 635 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 636* 637* Flag non-convergence. 638* 639 IF( J.GT.IOUT-IINFO ) THEN 640 NCNVRG = .TRUE. 641 IB = -JB 642 ELSE 643 IB = JB 644 END IF 645 DO 50 JE = IWORK( J ) + 1 + IWOFF, 646 $ IWORK( J+IN ) + IWOFF 647 W( JE ) = TMP1 648 IBLOCK( JE ) = IB 649 50 CONTINUE 650 60 CONTINUE 651* 652 M = M + IM 653 END IF 654 70 CONTINUE 655* 656* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 657* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 658* 659 IF( IRANGE.EQ.3 ) THEN 660 IM = 0 661 IDISCL = IL - 1 - NWL 662 IDISCU = NWU - IU 663* 664 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 665 DO 80 JE = 1, M 666 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 667 IDISCL = IDISCL - 1 668 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 669 IDISCU = IDISCU - 1 670 ELSE 671 IM = IM + 1 672 W( IM ) = W( JE ) 673 IBLOCK( IM ) = IBLOCK( JE ) 674 END IF 675 80 CONTINUE 676 M = IM 677 END IF 678 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 679* 680* Code to deal with effects of bad arithmetic: 681* Some low eigenvalues to be discarded are not in (WL,WLU], 682* or high eigenvalues to be discarded are not in (WUL,WU] 683* so just kill off the smallest IDISCL/largest IDISCU 684* eigenvalues, by simply finding the smallest/largest 685* eigenvalue(s). 686* 687* (If N(w) is monotone non-decreasing, this should never 688* happen.) 689* 690 IF( IDISCL.GT.0 ) THEN 691 WKILL = WU 692 DO 100 JDISC = 1, IDISCL 693 IW = 0 694 DO 90 JE = 1, M 695 IF( IBLOCK( JE ).NE.0 .AND. 696 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 697 IW = JE 698 WKILL = W( JE ) 699 END IF 700 90 CONTINUE 701 IBLOCK( IW ) = 0 702 100 CONTINUE 703 END IF 704 IF( IDISCU.GT.0 ) THEN 705* 706 WKILL = WL 707 DO 120 JDISC = 1, IDISCU 708 IW = 0 709 DO 110 JE = 1, M 710 IF( IBLOCK( JE ).NE.0 .AND. 711 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN 712 IW = JE 713 WKILL = W( JE ) 714 END IF 715 110 CONTINUE 716 IBLOCK( IW ) = 0 717 120 CONTINUE 718 END IF 719 IM = 0 720 DO 130 JE = 1, M 721 IF( IBLOCK( JE ).NE.0 ) THEN 722 IM = IM + 1 723 W( IM ) = W( JE ) 724 IBLOCK( IM ) = IBLOCK( JE ) 725 END IF 726 130 CONTINUE 727 M = IM 728 END IF 729 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 730 TOOFEW = .TRUE. 731 END IF 732 END IF 733* 734* If ORDER='B', do nothing -- the eigenvalues are already sorted 735* by block. 736* If ORDER='E', sort the eigenvalues from smallest to largest 737* 738 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN 739 DO 150 JE = 1, M - 1 740 IE = 0 741 TMP1 = W( JE ) 742 DO 140 J = JE + 1, M 743 IF( W( J ).LT.TMP1 ) THEN 744 IE = J 745 TMP1 = W( J ) 746 END IF 747 140 CONTINUE 748* 749 IF( IE.NE.0 ) THEN 750 ITMP1 = IBLOCK( IE ) 751 W( IE ) = W( JE ) 752 IBLOCK( IE ) = IBLOCK( JE ) 753 W( JE ) = TMP1 754 IBLOCK( JE ) = ITMP1 755 END IF 756 150 CONTINUE 757 END IF 758* 759 INFO = 0 760 IF( NCNVRG ) 761 $ INFO = INFO + 1 762 IF( TOOFEW ) 763 $ INFO = INFO + 2 764 RETURN 765* 766* End of DSTEBZ 767* 768 END 769