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