Annotation of rpl/lapack/lapack/dlarrd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
! 2: $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
! 3: $ M, W, WERR, WL, WU, IBLOCK, INDEXW,
! 4: $ WORK, IWORK, INFO )
! 5: *
! 6: * -- LAPACK auxiliary routine (version 3.2.1) --
! 7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 9: * -- April 2009 --
! 10: *
! 11: * .. Scalar Arguments ..
! 12: CHARACTER ORDER, RANGE
! 13: INTEGER IL, INFO, IU, M, N, NSPLIT
! 14: DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
! 15: * ..
! 16: * .. Array Arguments ..
! 17: INTEGER IBLOCK( * ), INDEXW( * ),
! 18: $ ISPLIT( * ), IWORK( * )
! 19: DOUBLE PRECISION D( * ), E( * ), E2( * ),
! 20: $ GERS( * ), W( * ), WERR( * ), WORK( * )
! 21: * ..
! 22: *
! 23: * Purpose
! 24: * =======
! 25: *
! 26: * DLARRD computes the eigenvalues of a symmetric tridiagonal
! 27: * matrix T to suitable accuracy. This is an auxiliary code to be
! 28: * called from DSTEMR.
! 29: * The user may ask for all eigenvalues, all eigenvalues
! 30: * in the half-open interval (VL, VU], or the IL-th through IU-th
! 31: * eigenvalues.
! 32: *
! 33: * To avoid overflow, the matrix must be scaled so that its
! 34: * largest element is no greater than overflow**(1/2) *
! 35: * underflow**(1/4) in absolute value, and for greatest
! 36: * accuracy, it should not be much smaller than that.
! 37: *
! 38: * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
! 39: * Matrix", Report CS41, Computer Science Dept., Stanford
! 40: * University, July 21, 1966.
! 41: *
! 42: * Arguments
! 43: * =========
! 44: *
! 45: * RANGE (input) CHARACTER
! 46: * = 'A': ("All") all eigenvalues will be found.
! 47: * = 'V': ("Value") all eigenvalues in the half-open interval
! 48: * (VL, VU] will be found.
! 49: * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
! 50: * entire matrix) will be found.
! 51: *
! 52: * ORDER (input) CHARACTER
! 53: * = 'B': ("By Block") the eigenvalues will be grouped by
! 54: * split-off block (see IBLOCK, ISPLIT) and
! 55: * ordered from smallest to largest within
! 56: * the block.
! 57: * = 'E': ("Entire matrix")
! 58: * the eigenvalues for the entire matrix
! 59: * will be ordered from smallest to
! 60: * largest.
! 61: *
! 62: * N (input) INTEGER
! 63: * The order of the tridiagonal matrix T. N >= 0.
! 64: *
! 65: * VL (input) DOUBLE PRECISION
! 66: * VU (input) DOUBLE PRECISION
! 67: * If RANGE='V', the lower and upper bounds of the interval to
! 68: * be searched for eigenvalues. Eigenvalues less than or equal
! 69: * to VL, or greater than VU, will not be returned. VL < VU.
! 70: * Not referenced if RANGE = 'A' or 'I'.
! 71: *
! 72: * IL (input) INTEGER
! 73: * IU (input) INTEGER
! 74: * If RANGE='I', the indices (in ascending order) of the
! 75: * smallest and largest eigenvalues to be returned.
! 76: * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
! 77: * Not referenced if RANGE = 'A' or 'V'.
! 78: *
! 79: * GERS (input) DOUBLE PRECISION array, dimension (2*N)
! 80: * The N Gerschgorin intervals (the i-th Gerschgorin interval
! 81: * is (GERS(2*i-1), GERS(2*i)).
! 82: *
! 83: * RELTOL (input) DOUBLE PRECISION
! 84: * The minimum relative width of an interval. When an interval
! 85: * is narrower than RELTOL times the larger (in
! 86: * magnitude) endpoint, then it is considered to be
! 87: * sufficiently small, i.e., converged. Note: this should
! 88: * always be at least radix*machine epsilon.
! 89: *
! 90: * D (input) DOUBLE PRECISION array, dimension (N)
! 91: * The n diagonal elements of the tridiagonal matrix T.
! 92: *
! 93: * E (input) DOUBLE PRECISION array, dimension (N-1)
! 94: * The (n-1) off-diagonal elements of the tridiagonal matrix T.
! 95: *
! 96: * E2 (input) DOUBLE PRECISION array, dimension (N-1)
! 97: * The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
! 98: *
! 99: * PIVMIN (input) DOUBLE PRECISION
! 100: * The minimum pivot allowed in the Sturm sequence for T.
! 101: *
! 102: * NSPLIT (input) INTEGER
! 103: * The number of diagonal blocks in the matrix T.
! 104: * 1 <= NSPLIT <= N.
! 105: *
! 106: * ISPLIT (input) INTEGER array, dimension (N)
! 107: * The splitting points, at which T breaks up into submatrices.
! 108: * The first submatrix consists of rows/columns 1 to ISPLIT(1),
! 109: * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
! 110: * etc., and the NSPLIT-th consists of rows/columns
! 111: * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
! 112: * (Only the first NSPLIT elements will actually be used, but
! 113: * since the user cannot know a priori what value NSPLIT will
! 114: * have, N words must be reserved for ISPLIT.)
! 115: *
! 116: * M (output) INTEGER
! 117: * The actual number of eigenvalues found. 0 <= M <= N.
! 118: * (See also the description of INFO=2,3.)
! 119: *
! 120: * W (output) DOUBLE PRECISION array, dimension (N)
! 121: * On exit, the first M elements of W will contain the
! 122: * eigenvalue approximations. DLARRD computes an interval
! 123: * I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
! 124: * approximation is given as the interval midpoint
! 125: * W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
! 126: * WERR(j) = abs( a_j - b_j)/2
! 127: *
! 128: * WERR (output) DOUBLE PRECISION array, dimension (N)
! 129: * The error bound on the corresponding eigenvalue approximation
! 130: * in W.
! 131: *
! 132: * WL (output) DOUBLE PRECISION
! 133: * WU (output) DOUBLE PRECISION
! 134: * The interval (WL, WU] contains all the wanted eigenvalues.
! 135: * If RANGE='V', then WL=VL and WU=VU.
! 136: * If RANGE='A', then WL and WU are the global Gerschgorin bounds
! 137: * on the spectrum.
! 138: * If RANGE='I', then WL and WU are computed by DLAEBZ from the
! 139: * index range specified.
! 140: *
! 141: * IBLOCK (output) INTEGER array, dimension (N)
! 142: * At each row/column j where E(j) is zero or small, the
! 143: * matrix T is considered to split into a block diagonal
! 144: * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
! 145: * block (from 1 to the number of blocks) the eigenvalue W(i)
! 146: * belongs. (DLARRD may use the remaining N-M elements as
! 147: * workspace.)
! 148: *
! 149: * INDEXW (output) INTEGER array, dimension (N)
! 150: * The indices of the eigenvalues within each block (submatrix);
! 151: * for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
! 152: * i-th eigenvalue W(i) is the j-th eigenvalue in block k.
! 153: *
! 154: * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
! 155: *
! 156: * IWORK (workspace) INTEGER array, dimension (3*N)
! 157: *
! 158: * INFO (output) INTEGER
! 159: * = 0: successful exit
! 160: * < 0: if INFO = -i, the i-th argument had an illegal value
! 161: * > 0: some or all of the eigenvalues failed to converge or
! 162: * were not computed:
! 163: * =1 or 3: Bisection failed to converge for some
! 164: * eigenvalues; these eigenvalues are flagged by a
! 165: * negative block number. The effect is that the
! 166: * eigenvalues may not be as accurate as the
! 167: * absolute and relative tolerances. This is
! 168: * generally caused by unexpectedly inaccurate
! 169: * arithmetic.
! 170: * =2 or 3: RANGE='I' only: Not all of the eigenvalues
! 171: * IL:IU were found.
! 172: * Effect: M < IU+1-IL
! 173: * Cause: non-monotonic arithmetic, causing the
! 174: * Sturm sequence to be non-monotonic.
! 175: * Cure: recalculate, using RANGE='A', and pick
! 176: * out eigenvalues IL:IU. In some cases,
! 177: * increasing the PARAMETER "FUDGE" may
! 178: * make things work.
! 179: * = 4: RANGE='I', and the Gershgorin interval
! 180: * initially used was too small. No eigenvalues
! 181: * were computed.
! 182: * Probable cause: your machine has sloppy
! 183: * floating-point arithmetic.
! 184: * Cure: Increase the PARAMETER "FUDGE",
! 185: * recompile, and try again.
! 186: *
! 187: * Internal Parameters
! 188: * ===================
! 189: *
! 190: * FUDGE DOUBLE PRECISION, default = 2
! 191: * A "fudge factor" to widen the Gershgorin intervals. Ideally,
! 192: * a value of 1 should work, but on machines with sloppy
! 193: * arithmetic, this needs to be larger. The default for
! 194: * publicly released versions should be large enough to handle
! 195: * the worst machine around. Note that this has no effect
! 196: * on accuracy of the solution.
! 197: *
! 198: * Based on contributions by
! 199: * W. Kahan, University of California, Berkeley, USA
! 200: * Beresford Parlett, University of California, Berkeley, USA
! 201: * Jim Demmel, University of California, Berkeley, USA
! 202: * Inderjit Dhillon, University of Texas, Austin, USA
! 203: * Osni Marques, LBNL/NERSC, USA
! 204: * Christof Voemel, University of California, Berkeley, USA
! 205: *
! 206: * =====================================================================
! 207: *
! 208: * .. Parameters ..
! 209: DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
! 210: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
! 211: $ TWO = 2.0D0, HALF = ONE/TWO,
! 212: $ FUDGE = TWO )
! 213: INTEGER ALLRNG, VALRNG, INDRNG
! 214: PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 )
! 215: * ..
! 216: * .. Local Scalars ..
! 217: LOGICAL NCNVRG, TOOFEW
! 218: INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
! 219: $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
! 220: $ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB,
! 221: $ NWL, NWU
! 222: DOUBLE PRECISION ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2,
! 223: $ TNORM, UFLOW, WKILL, WLU, WUL
! 224:
! 225: * ..
! 226: * .. Local Arrays ..
! 227: INTEGER IDUMMA( 1 )
! 228: * ..
! 229: * .. External Functions ..
! 230: LOGICAL LSAME
! 231: INTEGER ILAENV
! 232: DOUBLE PRECISION DLAMCH
! 233: EXTERNAL LSAME, ILAENV, DLAMCH
! 234: * ..
! 235: * .. External Subroutines ..
! 236: EXTERNAL DLAEBZ
! 237: * ..
! 238: * .. Intrinsic Functions ..
! 239: INTRINSIC ABS, INT, LOG, MAX, MIN
! 240: * ..
! 241: * .. Executable Statements ..
! 242: *
! 243: INFO = 0
! 244: *
! 245: * Decode RANGE
! 246: *
! 247: IF( LSAME( RANGE, 'A' ) ) THEN
! 248: IRANGE = ALLRNG
! 249: ELSE IF( LSAME( RANGE, 'V' ) ) THEN
! 250: IRANGE = VALRNG
! 251: ELSE IF( LSAME( RANGE, 'I' ) ) THEN
! 252: IRANGE = INDRNG
! 253: ELSE
! 254: IRANGE = 0
! 255: END IF
! 256: *
! 257: * Check for Errors
! 258: *
! 259: IF( IRANGE.LE.0 ) THEN
! 260: INFO = -1
! 261: ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN
! 262: INFO = -2
! 263: ELSE IF( N.LT.0 ) THEN
! 264: INFO = -3
! 265: ELSE IF( IRANGE.EQ.VALRNG ) THEN
! 266: IF( VL.GE.VU )
! 267: $ INFO = -5
! 268: ELSE IF( IRANGE.EQ.INDRNG .AND.
! 269: $ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN
! 270: INFO = -6
! 271: ELSE IF( IRANGE.EQ.INDRNG .AND.
! 272: $ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN
! 273: INFO = -7
! 274: END IF
! 275: *
! 276: IF( INFO.NE.0 ) THEN
! 277: RETURN
! 278: END IF
! 279:
! 280: * Initialize error flags
! 281: INFO = 0
! 282: NCNVRG = .FALSE.
! 283: TOOFEW = .FALSE.
! 284:
! 285: * Quick return if possible
! 286: M = 0
! 287: IF( N.EQ.0 ) RETURN
! 288:
! 289: * Simplification:
! 290: IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1
! 291:
! 292: * Get machine constants
! 293: EPS = DLAMCH( 'P' )
! 294: UFLOW = DLAMCH( 'U' )
! 295:
! 296:
! 297: * Special Case when N=1
! 298: * Treat case of 1x1 matrix for quick return
! 299: IF( N.EQ.1 ) THEN
! 300: IF( (IRANGE.EQ.ALLRNG).OR.
! 301: $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
! 302: $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
! 303: M = 1
! 304: W(1) = D(1)
! 305: * The computation error of the eigenvalue is zero
! 306: WERR(1) = ZERO
! 307: IBLOCK( 1 ) = 1
! 308: INDEXW( 1 ) = 1
! 309: ENDIF
! 310: RETURN
! 311: END IF
! 312:
! 313: * NB is the minimum vector length for vector bisection, or 0
! 314: * if only scalar is to be done.
! 315: NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
! 316: IF( NB.LE.1 ) NB = 0
! 317:
! 318: * Find global spectral radius
! 319: GL = D(1)
! 320: GU = D(1)
! 321: DO 5 I = 1,N
! 322: GL = MIN( GL, GERS( 2*I - 1))
! 323: GU = MAX( GU, GERS(2*I) )
! 324: 5 CONTINUE
! 325: * Compute global Gerschgorin bounds and spectral diameter
! 326: TNORM = MAX( ABS( GL ), ABS( GU ) )
! 327: GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
! 328: GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
! 329: * [JAN/28/2009] remove the line below since SPDIAM variable not use
! 330: * SPDIAM = GU - GL
! 331: * Input arguments for DLAEBZ:
! 332: * The relative tolerance. An interval (a,b] lies within
! 333: * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
! 334: RTOLI = RELTOL
! 335: * Set the absolute tolerance for interval convergence to zero to force
! 336: * interval convergence based on relative size of the interval.
! 337: * This is dangerous because intervals might not converge when RELTOL is
! 338: * small. But at least a very small number should be selected so that for
! 339: * strongly graded matrices, the code can get relatively accurate
! 340: * eigenvalues.
! 341: ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN
! 342:
! 343: IF( IRANGE.EQ.INDRNG ) THEN
! 344:
! 345: * RANGE='I': Compute an interval containing eigenvalues
! 346: * IL through IU. The initial interval [GL,GU] from the global
! 347: * Gerschgorin bounds GL and GU is refined by DLAEBZ.
! 348: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
! 349: $ LOG( TWO ) ) + 2
! 350: WORK( N+1 ) = GL
! 351: WORK( N+2 ) = GL
! 352: WORK( N+3 ) = GU
! 353: WORK( N+4 ) = GU
! 354: WORK( N+5 ) = GL
! 355: WORK( N+6 ) = GU
! 356: IWORK( 1 ) = -1
! 357: IWORK( 2 ) = -1
! 358: IWORK( 3 ) = N + 1
! 359: IWORK( 4 ) = N + 1
! 360: IWORK( 5 ) = IL - 1
! 361: IWORK( 6 ) = IU
! 362: *
! 363: CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN,
! 364: $ D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
! 365: $ IWORK, W, IBLOCK, IINFO )
! 366: IF( IINFO .NE. 0 ) THEN
! 367: INFO = IINFO
! 368: RETURN
! 369: END IF
! 370: * On exit, output intervals may not be ordered by ascending negcount
! 371: IF( IWORK( 6 ).EQ.IU ) THEN
! 372: WL = WORK( N+1 )
! 373: WLU = WORK( N+3 )
! 374: NWL = IWORK( 1 )
! 375: WU = WORK( N+4 )
! 376: WUL = WORK( N+2 )
! 377: NWU = IWORK( 4 )
! 378: ELSE
! 379: WL = WORK( N+2 )
! 380: WLU = WORK( N+4 )
! 381: NWL = IWORK( 2 )
! 382: WU = WORK( N+3 )
! 383: WUL = WORK( N+1 )
! 384: NWU = IWORK( 3 )
! 385: END IF
! 386: * On exit, the interval [WL, WLU] contains a value with negcount NWL,
! 387: * and [WUL, WU] contains a value with negcount NWU.
! 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:
! 393: ELSEIF( IRANGE.EQ.VALRNG ) THEN
! 394: WL = VL
! 395: WU = VU
! 396:
! 397: ELSEIF( IRANGE.EQ.ALLRNG ) THEN
! 398: WL = GL
! 399: WU = GU
! 400: ENDIF
! 401:
! 402:
! 403:
! 404: * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
! 405: * NWL accumulates the number of eigenvalues .le. WL,
! 406: * NWU accumulates the number of eigenvalues .le. WU
! 407: M = 0
! 408: IEND = 0
! 409: INFO = 0
! 410: NWL = 0
! 411: NWU = 0
! 412: *
! 413: DO 70 JBLK = 1, NSPLIT
! 414: IOFF = IEND
! 415: IBEGIN = IOFF + 1
! 416: IEND = ISPLIT( JBLK )
! 417: IN = IEND - IOFF
! 418: *
! 419: IF( IN.EQ.1 ) THEN
! 420: * 1x1 block
! 421: IF( WL.GE.D( IBEGIN )-PIVMIN )
! 422: $ NWL = NWL + 1
! 423: IF( WU.GE.D( IBEGIN )-PIVMIN )
! 424: $ NWU = NWU + 1
! 425: IF( IRANGE.EQ.ALLRNG .OR.
! 426: $ ( WL.LT.D( IBEGIN )-PIVMIN
! 427: $ .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN
! 428: M = M + 1
! 429: W( M ) = D( IBEGIN )
! 430: WERR(M) = ZERO
! 431: * The gap for a single block doesn't matter for the later
! 432: * algorithm and is assigned an arbitrary large value
! 433: IBLOCK( M ) = JBLK
! 434: INDEXW( M ) = 1
! 435: END IF
! 436:
! 437: * Disabled 2x2 case because of a failure on the following matrix
! 438: * RANGE = 'I', IL = IU = 4
! 439: * Original Tridiagonal, d = [
! 440: * -0.150102010615740E+00
! 441: * -0.849897989384260E+00
! 442: * -0.128208148052635E-15
! 443: * 0.128257718286320E-15
! 444: * ];
! 445: * e = [
! 446: * -0.357171383266986E+00
! 447: * -0.180411241501588E-15
! 448: * -0.175152352710251E-15
! 449: * ];
! 450: *
! 451: * ELSE IF( IN.EQ.2 ) THEN
! 452: ** 2x2 block
! 453: * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
! 454: * TMP1 = HALF*(D(IBEGIN)+D(IEND))
! 455: * L1 = TMP1 - DISC
! 456: * IF( WL.GE. L1-PIVMIN )
! 457: * $ NWL = NWL + 1
! 458: * IF( WU.GE. L1-PIVMIN )
! 459: * $ NWU = NWU + 1
! 460: * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
! 461: * $ L1-PIVMIN ) ) THEN
! 462: * M = M + 1
! 463: * W( M ) = L1
! 464: ** The uncertainty of eigenvalues of a 2x2 matrix is very small
! 465: * WERR( M ) = EPS * ABS( W( M ) ) * TWO
! 466: * IBLOCK( M ) = JBLK
! 467: * INDEXW( M ) = 1
! 468: * ENDIF
! 469: * L2 = TMP1 + DISC
! 470: * IF( WL.GE. L2-PIVMIN )
! 471: * $ NWL = NWL + 1
! 472: * IF( WU.GE. L2-PIVMIN )
! 473: * $ NWU = NWU + 1
! 474: * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
! 475: * $ L2-PIVMIN ) ) THEN
! 476: * M = M + 1
! 477: * W( M ) = L2
! 478: ** The uncertainty of eigenvalues of a 2x2 matrix is very small
! 479: * WERR( M ) = EPS * ABS( W( M ) ) * TWO
! 480: * IBLOCK( M ) = JBLK
! 481: * INDEXW( M ) = 2
! 482: * ENDIF
! 483: ELSE
! 484: * General Case - block of size IN >= 2
! 485: * Compute local Gerschgorin interval and use it as the initial
! 486: * interval for DLAEBZ
! 487: GU = D( IBEGIN )
! 488: GL = D( IBEGIN )
! 489: TMP1 = ZERO
! 490:
! 491: DO 40 J = IBEGIN, IEND
! 492: GL = MIN( GL, GERS( 2*J - 1))
! 493: GU = MAX( GU, GERS(2*J) )
! 494: 40 CONTINUE
! 495: * [JAN/28/2009]
! 496: * change SPDIAM by TNORM in lines 2 and 3 thereafter
! 497: * line 1: remove computation of SPDIAM (not useful anymore)
! 498: * SPDIAM = GU - GL
! 499: * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
! 500: * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
! 501: GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN
! 502: GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN
! 503: *
! 504: IF( IRANGE.GT.1 ) THEN
! 505: IF( GU.LT.WL ) THEN
! 506: * the local block contains none of the wanted eigenvalues
! 507: NWL = NWL + IN
! 508: NWU = NWU + IN
! 509: GO TO 70
! 510: END IF
! 511: * refine search interval if possible, only range (WL,WU] matters
! 512: GL = MAX( GL, WL )
! 513: GU = MIN( GU, WU )
! 514: IF( GL.GE.GU )
! 515: $ GO TO 70
! 516: END IF
! 517:
! 518: * Find negcount of initial interval boundaries GL and GU
! 519: WORK( N+1 ) = GL
! 520: WORK( N+IN+1 ) = GU
! 521: CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
! 522: $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
! 523: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
! 524: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
! 525: IF( IINFO .NE. 0 ) THEN
! 526: INFO = IINFO
! 527: RETURN
! 528: END IF
! 529: *
! 530: NWL = NWL + IWORK( 1 )
! 531: NWU = NWU + IWORK( IN+1 )
! 532: IWOFF = M - IWORK( 1 )
! 533:
! 534: * Compute Eigenvalues
! 535: ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
! 536: $ LOG( TWO ) ) + 2
! 537: CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
! 538: $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ),
! 539: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
! 540: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
! 541: IF( IINFO .NE. 0 ) THEN
! 542: INFO = IINFO
! 543: RETURN
! 544: END IF
! 545: *
! 546: * Copy eigenvalues into W and IBLOCK
! 547: * Use -JBLK for block number for unconverged eigenvalues.
! 548: * Loop over the number of output intervals from DLAEBZ
! 549: DO 60 J = 1, IOUT
! 550: * eigenvalue approximation is middle point of interval
! 551: TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
! 552: * semi length of error interval
! 553: TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) )
! 554: IF( J.GT.IOUT-IINFO ) THEN
! 555: * Flag non-convergence.
! 556: NCNVRG = .TRUE.
! 557: IB = -JBLK
! 558: ELSE
! 559: IB = JBLK
! 560: END IF
! 561: DO 50 JE = IWORK( J ) + 1 + IWOFF,
! 562: $ IWORK( J+IN ) + IWOFF
! 563: W( JE ) = TMP1
! 564: WERR( JE ) = TMP2
! 565: INDEXW( JE ) = JE - IWOFF
! 566: IBLOCK( JE ) = IB
! 567: 50 CONTINUE
! 568: 60 CONTINUE
! 569: *
! 570: M = M + IM
! 571: END IF
! 572: 70 CONTINUE
! 573:
! 574: * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
! 575: * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
! 576: IF( IRANGE.EQ.INDRNG ) THEN
! 577: IDISCL = IL - 1 - NWL
! 578: IDISCU = NWU - IU
! 579: *
! 580: IF( IDISCL.GT.0 ) THEN
! 581: IM = 0
! 582: DO 80 JE = 1, M
! 583: * Remove some of the smallest eigenvalues from the left so that
! 584: * at the end IDISCL =0. Move all eigenvalues up to the left.
! 585: IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
! 586: IDISCL = IDISCL - 1
! 587: ELSE
! 588: IM = IM + 1
! 589: W( IM ) = W( JE )
! 590: WERR( IM ) = WERR( JE )
! 591: INDEXW( IM ) = INDEXW( JE )
! 592: IBLOCK( IM ) = IBLOCK( JE )
! 593: END IF
! 594: 80 CONTINUE
! 595: M = IM
! 596: END IF
! 597: IF( IDISCU.GT.0 ) THEN
! 598: * Remove some of the largest eigenvalues from the right so that
! 599: * at the end IDISCU =0. Move all eigenvalues up to the left.
! 600: IM=M+1
! 601: DO 81 JE = M, 1, -1
! 602: IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
! 603: IDISCU = IDISCU - 1
! 604: ELSE
! 605: IM = IM - 1
! 606: W( IM ) = W( JE )
! 607: WERR( IM ) = WERR( JE )
! 608: INDEXW( IM ) = INDEXW( JE )
! 609: IBLOCK( IM ) = IBLOCK( JE )
! 610: END IF
! 611: 81 CONTINUE
! 612: JEE = 0
! 613: DO 82 JE = IM, M
! 614: JEE = JEE + 1
! 615: W( JEE ) = W( JE )
! 616: WERR( JEE ) = WERR( JE )
! 617: INDEXW( JEE ) = INDEXW( JE )
! 618: IBLOCK( JEE ) = IBLOCK( JE )
! 619: 82 CONTINUE
! 620: M = M-IM+1
! 621: END IF
! 622:
! 623: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
! 624: * Code to deal with effects of bad arithmetic. (If N(w) is
! 625: * monotone non-decreasing, this should never happen.)
! 626: * Some low eigenvalues to be discarded are not in (WL,WLU],
! 627: * or high eigenvalues to be discarded are not in (WUL,WU]
! 628: * so just kill off the smallest IDISCL/largest IDISCU
! 629: * eigenvalues, by marking the corresponding IBLOCK = 0
! 630: IF( IDISCL.GT.0 ) THEN
! 631: WKILL = WU
! 632: DO 100 JDISC = 1, IDISCL
! 633: IW = 0
! 634: DO 90 JE = 1, M
! 635: IF( IBLOCK( JE ).NE.0 .AND.
! 636: $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
! 637: IW = JE
! 638: WKILL = W( JE )
! 639: END IF
! 640: 90 CONTINUE
! 641: IBLOCK( IW ) = 0
! 642: 100 CONTINUE
! 643: END IF
! 644: IF( IDISCU.GT.0 ) THEN
! 645: WKILL = WL
! 646: DO 120 JDISC = 1, IDISCU
! 647: IW = 0
! 648: DO 110 JE = 1, M
! 649: IF( IBLOCK( JE ).NE.0 .AND.
! 650: $ ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN
! 651: IW = JE
! 652: WKILL = W( JE )
! 653: END IF
! 654: 110 CONTINUE
! 655: IBLOCK( IW ) = 0
! 656: 120 CONTINUE
! 657: END IF
! 658: * Now erase all eigenvalues with IBLOCK set to zero
! 659: IM = 0
! 660: DO 130 JE = 1, M
! 661: IF( IBLOCK( JE ).NE.0 ) THEN
! 662: IM = IM + 1
! 663: W( IM ) = W( JE )
! 664: WERR( IM ) = WERR( JE )
! 665: INDEXW( IM ) = INDEXW( JE )
! 666: IBLOCK( IM ) = IBLOCK( JE )
! 667: END IF
! 668: 130 CONTINUE
! 669: M = IM
! 670: END IF
! 671: IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
! 672: TOOFEW = .TRUE.
! 673: END IF
! 674: END IF
! 675: *
! 676: IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR.
! 677: $ ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN
! 678: TOOFEW = .TRUE.
! 679: END IF
! 680:
! 681: * If ORDER='B', do nothing the eigenvalues are already sorted by
! 682: * block.
! 683: * If ORDER='E', sort the eigenvalues from smallest to largest
! 684:
! 685: IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN
! 686: DO 150 JE = 1, M - 1
! 687: IE = 0
! 688: TMP1 = W( JE )
! 689: DO 140 J = JE + 1, M
! 690: IF( W( J ).LT.TMP1 ) THEN
! 691: IE = J
! 692: TMP1 = W( J )
! 693: END IF
! 694: 140 CONTINUE
! 695: IF( IE.NE.0 ) THEN
! 696: TMP2 = WERR( IE )
! 697: ITMP1 = IBLOCK( IE )
! 698: ITMP2 = INDEXW( IE )
! 699: W( IE ) = W( JE )
! 700: WERR( IE ) = WERR( JE )
! 701: IBLOCK( IE ) = IBLOCK( JE )
! 702: INDEXW( IE ) = INDEXW( JE )
! 703: W( JE ) = TMP1
! 704: WERR( JE ) = TMP2
! 705: IBLOCK( JE ) = ITMP1
! 706: INDEXW( JE ) = ITMP2
! 707: END IF
! 708: 150 CONTINUE
! 709: END IF
! 710: *
! 711: INFO = 0
! 712: IF( NCNVRG )
! 713: $ INFO = INFO + 1
! 714: IF( TOOFEW )
! 715: $ INFO = INFO + 2
! 716: RETURN
! 717: *
! 718: * End of DLARRD
! 719: *
! 720: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack