Annotation of rpl/lapack/lapack/dstebz.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dstebz.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack