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