Annotation of rpl/lapack/lapack/dlarre.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
! 2: $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
! 3: $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
! 4: $ WORK, IWORK, INFO )
! 5: IMPLICIT NONE
! 6: *
! 7: * -- LAPACK auxiliary routine (version 3.2) --
! 8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 10: * November 2006
! 11: *
! 12: * .. Scalar Arguments ..
! 13: CHARACTER RANGE
! 14: INTEGER IL, INFO, IU, M, N, NSPLIT
! 15: DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
! 16: * ..
! 17: * .. Array Arguments ..
! 18: INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
! 19: $ INDEXW( * )
! 20: DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
! 21: $ W( * ),WERR( * ), WGAP( * ), WORK( * )
! 22: * ..
! 23: *
! 24: * Purpose
! 25: * =======
! 26: *
! 27: * To find the desired eigenvalues of a given real symmetric
! 28: * tridiagonal matrix T, DLARRE sets any "small" off-diagonal
! 29: * elements to zero, and for each unreduced block T_i, it finds
! 30: * (a) a suitable shift at one end of the block's spectrum,
! 31: * (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
! 32: * (c) eigenvalues of each L_i D_i L_i^T.
! 33: * The representations and eigenvalues found are then used by
! 34: * DSTEMR to compute the eigenvectors of T.
! 35: * The accuracy varies depending on whether bisection is used to
! 36: * find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
! 37: * conpute all and then discard any unwanted one.
! 38: * As an added benefit, DLARRE also outputs the n
! 39: * Gerschgorin intervals for the matrices L_i D_i L_i^T.
! 40: *
! 41: * Arguments
! 42: * =========
! 43: *
! 44: * RANGE (input) CHARACTER
! 45: * = 'A': ("All") all eigenvalues will be found.
! 46: * = 'V': ("Value") all eigenvalues in the half-open interval
! 47: * (VL, VU] will be found.
! 48: * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
! 49: * entire matrix) will be found.
! 50: *
! 51: * N (input) INTEGER
! 52: * The order of the matrix. N > 0.
! 53: *
! 54: * VL (input/output) DOUBLE PRECISION
! 55: * VU (input/output) DOUBLE PRECISION
! 56: * If RANGE='V', the lower and upper bounds for the eigenvalues.
! 57: * Eigenvalues less than or equal to VL, or greater than VU,
! 58: * will not be returned. VL < VU.
! 59: * If RANGE='I' or ='A', DLARRE computes bounds on the desired
! 60: * part of the spectrum.
! 61: *
! 62: * IL (input) INTEGER
! 63: * IU (input) INTEGER
! 64: * If RANGE='I', the indices (in ascending order) of the
! 65: * smallest and largest eigenvalues to be returned.
! 66: * 1 <= IL <= IU <= N.
! 67: *
! 68: * D (input/output) DOUBLE PRECISION array, dimension (N)
! 69: * On entry, the N diagonal elements of the tridiagonal
! 70: * matrix T.
! 71: * On exit, the N diagonal elements of the diagonal
! 72: * matrices D_i.
! 73: *
! 74: * E (input/output) DOUBLE PRECISION array, dimension (N)
! 75: * On entry, the first (N-1) entries contain the subdiagonal
! 76: * elements of the tridiagonal matrix T; E(N) need not be set.
! 77: * On exit, E contains the subdiagonal elements of the unit
! 78: * bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
! 79: * 1 <= I <= NSPLIT, contain the base points sigma_i on output.
! 80: *
! 81: * E2 (input/output) DOUBLE PRECISION array, dimension (N)
! 82: * On entry, the first (N-1) entries contain the SQUARES of the
! 83: * subdiagonal elements of the tridiagonal matrix T;
! 84: * E2(N) need not be set.
! 85: * On exit, the entries E2( ISPLIT( I ) ),
! 86: * 1 <= I <= NSPLIT, have been set to zero
! 87: *
! 88: * RTOL1 (input) DOUBLE PRECISION
! 89: * RTOL2 (input) DOUBLE PRECISION
! 90: * Parameters for bisection.
! 91: * An interval [LEFT,RIGHT] has converged if
! 92: * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
! 93: *
! 94: * SPLTOL (input) DOUBLE PRECISION
! 95: * The threshold for splitting.
! 96: *
! 97: * NSPLIT (output) INTEGER
! 98: * The number of blocks T splits into. 1 <= NSPLIT <= N.
! 99: *
! 100: * ISPLIT (output) INTEGER array, dimension (N)
! 101: * The splitting points, at which T breaks up into blocks.
! 102: * The first block consists of rows/columns 1 to ISPLIT(1),
! 103: * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
! 104: * etc., and the NSPLIT-th consists of rows/columns
! 105: * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
! 106: *
! 107: * M (output) INTEGER
! 108: * The total number of eigenvalues (of all L_i D_i L_i^T)
! 109: * found.
! 110: *
! 111: * W (output) DOUBLE PRECISION array, dimension (N)
! 112: * The first M elements contain the eigenvalues. The
! 113: * eigenvalues of each of the blocks, L_i D_i L_i^T, are
! 114: * sorted in ascending order ( DLARRE may use the
! 115: * remaining N-M elements as workspace).
! 116: *
! 117: * WERR (output) DOUBLE PRECISION array, dimension (N)
! 118: * The error bound on the corresponding eigenvalue in W.
! 119: *
! 120: * WGAP (output) DOUBLE PRECISION array, dimension (N)
! 121: * The separation from the right neighbor eigenvalue in W.
! 122: * The gap is only with respect to the eigenvalues of the same block
! 123: * as each block has its own representation tree.
! 124: * Exception: at the right end of a block we store the left gap
! 125: *
! 126: * IBLOCK (output) INTEGER array, dimension (N)
! 127: * The indices of the blocks (submatrices) associated with the
! 128: * corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
! 129: * W(i) belongs to the first block from the top, =2 if W(i)
! 130: * belongs to the second block, etc.
! 131: *
! 132: * INDEXW (output) INTEGER array, dimension (N)
! 133: * The indices of the eigenvalues within each block (submatrix);
! 134: * for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
! 135: * i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
! 136: *
! 137: * GERS (output) DOUBLE PRECISION array, dimension (2*N)
! 138: * The N Gerschgorin intervals (the i-th Gerschgorin interval
! 139: * is (GERS(2*i-1), GERS(2*i)).
! 140: *
! 141: * PIVMIN (output) DOUBLE PRECISION
! 142: * The minimum pivot in the Sturm sequence for T.
! 143: *
! 144: * WORK (workspace) DOUBLE PRECISION array, dimension (6*N)
! 145: * Workspace.
! 146: *
! 147: * IWORK (workspace) INTEGER array, dimension (5*N)
! 148: * Workspace.
! 149: *
! 150: * INFO (output) INTEGER
! 151: * = 0: successful exit
! 152: * > 0: A problem occured in DLARRE.
! 153: * < 0: One of the called subroutines signaled an internal problem.
! 154: * Needs inspection of the corresponding parameter IINFO
! 155: * for further information.
! 156: *
! 157: * =-1: Problem in DLARRD.
! 158: * = 2: No base representation could be found in MAXTRY iterations.
! 159: * Increasing MAXTRY and recompilation might be a remedy.
! 160: * =-3: Problem in DLARRB when computing the refined root
! 161: * representation for DLASQ2.
! 162: * =-4: Problem in DLARRB when preforming bisection on the
! 163: * desired part of the spectrum.
! 164: * =-5: Problem in DLASQ2.
! 165: * =-6: Problem in DLASQ2.
! 166: *
! 167: * Further Details
! 168: * The base representations are required to suffer very little
! 169: * element growth and consequently define all their eigenvalues to
! 170: * high relative accuracy.
! 171: * ===============
! 172: *
! 173: * Based on contributions by
! 174: * Beresford Parlett, University of California, Berkeley, USA
! 175: * Jim Demmel, University of California, Berkeley, USA
! 176: * Inderjit Dhillon, University of Texas, Austin, USA
! 177: * Osni Marques, LBNL/NERSC, USA
! 178: * Christof Voemel, University of California, Berkeley, USA
! 179: *
! 180: * =====================================================================
! 181: *
! 182: * .. Parameters ..
! 183: DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
! 184: $ MAXGROWTH, ONE, PERT, TWO, ZERO
! 185: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
! 186: $ TWO = 2.0D0, FOUR=4.0D0,
! 187: $ HNDRD = 100.0D0,
! 188: $ PERT = 8.0D0,
! 189: $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
! 190: $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 )
! 191: INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
! 192: PARAMETER ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
! 193: $ VALRNG = 3 )
! 194: * ..
! 195: * .. Local Scalars ..
! 196: LOGICAL FORCEB, NOREP, USEDQD
! 197: INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
! 198: $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
! 199: $ WBEGIN, WEND
! 200: DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
! 201: $ EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
! 202: $ RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
! 203: $ TAU, TMP, TMP1
! 204:
! 205:
! 206: * ..
! 207: * .. Local Arrays ..
! 208: INTEGER ISEED( 4 )
! 209: * ..
! 210: * .. External Functions ..
! 211: LOGICAL LSAME
! 212: DOUBLE PRECISION DLAMCH
! 213: EXTERNAL DLAMCH, LSAME
! 214:
! 215: * ..
! 216: * .. External Subroutines ..
! 217: EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD,
! 218: $ DLASQ2
! 219: * ..
! 220: * .. Intrinsic Functions ..
! 221: INTRINSIC ABS, MAX, MIN
! 222:
! 223: * ..
! 224: * .. Executable Statements ..
! 225: *
! 226:
! 227: INFO = 0
! 228:
! 229: *
! 230: * Decode RANGE
! 231: *
! 232: IF( LSAME( RANGE, 'A' ) ) THEN
! 233: IRANGE = ALLRNG
! 234: ELSE IF( LSAME( RANGE, 'V' ) ) THEN
! 235: IRANGE = VALRNG
! 236: ELSE IF( LSAME( RANGE, 'I' ) ) THEN
! 237: IRANGE = INDRNG
! 238: END IF
! 239:
! 240: M = 0
! 241:
! 242: * Get machine constants
! 243: SAFMIN = DLAMCH( 'S' )
! 244: EPS = DLAMCH( 'P' )
! 245:
! 246: * Set parameters
! 247: RTL = SQRT(EPS)
! 248: BSRTOL = SQRT(EPS)
! 249:
! 250: * Treat case of 1x1 matrix for quick return
! 251: IF( N.EQ.1 ) THEN
! 252: IF( (IRANGE.EQ.ALLRNG).OR.
! 253: $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
! 254: $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
! 255: M = 1
! 256: W(1) = D(1)
! 257: * The computation error of the eigenvalue is zero
! 258: WERR(1) = ZERO
! 259: WGAP(1) = ZERO
! 260: IBLOCK( 1 ) = 1
! 261: INDEXW( 1 ) = 1
! 262: GERS(1) = D( 1 )
! 263: GERS(2) = D( 1 )
! 264: ENDIF
! 265: * store the shift for the initial RRR, which is zero in this case
! 266: E(1) = ZERO
! 267: RETURN
! 268: END IF
! 269:
! 270: * General case: tridiagonal matrix of order > 1
! 271: *
! 272: * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
! 273: * Compute maximum off-diagonal entry and pivmin.
! 274: GL = D(1)
! 275: GU = D(1)
! 276: EOLD = ZERO
! 277: EMAX = ZERO
! 278: E(N) = ZERO
! 279: DO 5 I = 1,N
! 280: WERR(I) = ZERO
! 281: WGAP(I) = ZERO
! 282: EABS = ABS( E(I) )
! 283: IF( EABS .GE. EMAX ) THEN
! 284: EMAX = EABS
! 285: END IF
! 286: TMP1 = EABS + EOLD
! 287: GERS( 2*I-1) = D(I) - TMP1
! 288: GL = MIN( GL, GERS( 2*I - 1))
! 289: GERS( 2*I ) = D(I) + TMP1
! 290: GU = MAX( GU, GERS(2*I) )
! 291: EOLD = EABS
! 292: 5 CONTINUE
! 293: * The minimum pivot allowed in the Sturm sequence for T
! 294: PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
! 295: * Compute spectral diameter. The Gerschgorin bounds give an
! 296: * estimate that is wrong by at most a factor of SQRT(2)
! 297: SPDIAM = GU - GL
! 298:
! 299: * Compute splitting points
! 300: CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM,
! 301: $ NSPLIT, ISPLIT, IINFO )
! 302:
! 303: * Can force use of bisection instead of faster DQDS.
! 304: * Option left in the code for future multisection work.
! 305: FORCEB = .FALSE.
! 306:
! 307: * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
! 308: * explicitly wants bisection.
! 309: USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
! 310:
! 311: IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
! 312: * Set interval [VL,VU] that contains all eigenvalues
! 313: VL = GL
! 314: VU = GU
! 315: ELSE
! 316: * We call DLARRD to find crude approximations to the eigenvalues
! 317: * in the desired range. In case IRANGE = INDRNG, we also obtain the
! 318: * interval (VL,VU] that contains all the wanted eigenvalues.
! 319: * An interval [LEFT,RIGHT] has converged if
! 320: * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
! 321: * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
! 322: CALL DLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
! 323: $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
! 324: $ MM, W, WERR, VL, VU, IBLOCK, INDEXW,
! 325: $ WORK, IWORK, IINFO )
! 326: IF( IINFO.NE.0 ) THEN
! 327: INFO = -1
! 328: RETURN
! 329: ENDIF
! 330: * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
! 331: DO 14 I = MM+1,N
! 332: W( I ) = ZERO
! 333: WERR( I ) = ZERO
! 334: IBLOCK( I ) = 0
! 335: INDEXW( I ) = 0
! 336: 14 CONTINUE
! 337: END IF
! 338:
! 339:
! 340: ***
! 341: * Loop over unreduced blocks
! 342: IBEGIN = 1
! 343: WBEGIN = 1
! 344: DO 170 JBLK = 1, NSPLIT
! 345: IEND = ISPLIT( JBLK )
! 346: IN = IEND - IBEGIN + 1
! 347:
! 348: * 1 X 1 block
! 349: IF( IN.EQ.1 ) THEN
! 350: IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
! 351: $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
! 352: $ .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
! 353: $ ) THEN
! 354: M = M + 1
! 355: W( M ) = D( IBEGIN )
! 356: WERR(M) = ZERO
! 357: * The gap for a single block doesn't matter for the later
! 358: * algorithm and is assigned an arbitrary large value
! 359: WGAP(M) = ZERO
! 360: IBLOCK( M ) = JBLK
! 361: INDEXW( M ) = 1
! 362: WBEGIN = WBEGIN + 1
! 363: ENDIF
! 364: * E( IEND ) holds the shift for the initial RRR
! 365: E( IEND ) = ZERO
! 366: IBEGIN = IEND + 1
! 367: GO TO 170
! 368: END IF
! 369: *
! 370: * Blocks of size larger than 1x1
! 371: *
! 372: * E( IEND ) will hold the shift for the initial RRR, for now set it =0
! 373: E( IEND ) = ZERO
! 374: *
! 375: * Find local outer bounds GL,GU for the block
! 376: GL = D(IBEGIN)
! 377: GU = D(IBEGIN)
! 378: DO 15 I = IBEGIN , IEND
! 379: GL = MIN( GERS( 2*I-1 ), GL )
! 380: GU = MAX( GERS( 2*I ), GU )
! 381: 15 CONTINUE
! 382: SPDIAM = GU - GL
! 383:
! 384: IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
! 385: * Count the number of eigenvalues in the current block.
! 386: MB = 0
! 387: DO 20 I = WBEGIN,MM
! 388: IF( IBLOCK(I).EQ.JBLK ) THEN
! 389: MB = MB+1
! 390: ELSE
! 391: GOTO 21
! 392: ENDIF
! 393: 20 CONTINUE
! 394: 21 CONTINUE
! 395:
! 396: IF( MB.EQ.0) THEN
! 397: * No eigenvalue in the current block lies in the desired range
! 398: * E( IEND ) holds the shift for the initial RRR
! 399: E( IEND ) = ZERO
! 400: IBEGIN = IEND + 1
! 401: GO TO 170
! 402: ELSE
! 403:
! 404: * Decide whether dqds or bisection is more efficient
! 405: USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
! 406: WEND = WBEGIN + MB - 1
! 407: * Calculate gaps for the current block
! 408: * In later stages, when representations for individual
! 409: * eigenvalues are different, we use SIGMA = E( IEND ).
! 410: SIGMA = ZERO
! 411: DO 30 I = WBEGIN, WEND - 1
! 412: WGAP( I ) = MAX( ZERO,
! 413: $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
! 414: 30 CONTINUE
! 415: WGAP( WEND ) = MAX( ZERO,
! 416: $ VU - SIGMA - (W( WEND )+WERR( WEND )))
! 417: * Find local index of the first and last desired evalue.
! 418: INDL = INDEXW(WBEGIN)
! 419: INDU = INDEXW( WEND )
! 420: ENDIF
! 421: ENDIF
! 422: IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
! 423: * Case of DQDS
! 424: * Find approximations to the extremal eigenvalues of the block
! 425: CALL DLARRK( IN, 1, GL, GU, D(IBEGIN),
! 426: $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
! 427: IF( IINFO.NE.0 ) THEN
! 428: INFO = -1
! 429: RETURN
! 430: ENDIF
! 431: ISLEFT = MAX(GL, TMP - TMP1
! 432: $ - HNDRD * EPS* ABS(TMP - TMP1))
! 433:
! 434: CALL DLARRK( IN, IN, GL, GU, D(IBEGIN),
! 435: $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
! 436: IF( IINFO.NE.0 ) THEN
! 437: INFO = -1
! 438: RETURN
! 439: ENDIF
! 440: ISRGHT = MIN(GU, TMP + TMP1
! 441: $ + HNDRD * EPS * ABS(TMP + TMP1))
! 442: * Improve the estimate of the spectral diameter
! 443: SPDIAM = ISRGHT - ISLEFT
! 444: ELSE
! 445: * Case of bisection
! 446: * Find approximations to the wanted extremal eigenvalues
! 447: ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
! 448: $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
! 449: ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
! 450: $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
! 451: ENDIF
! 452:
! 453:
! 454: * Decide whether the base representation for the current block
! 455: * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
! 456: * should be on the left or the right end of the current block.
! 457: * The strategy is to shift to the end which is "more populated"
! 458: * Furthermore, decide whether to use DQDS for the computation of
! 459: * the eigenvalue approximations at the end of DLARRE or bisection.
! 460: * dqds is chosen if all eigenvalues are desired or the number of
! 461: * eigenvalues to be computed is large compared to the blocksize.
! 462: IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
! 463: * If all the eigenvalues have to be computed, we use dqd
! 464: USEDQD = .TRUE.
! 465: * INDL is the local index of the first eigenvalue to compute
! 466: INDL = 1
! 467: INDU = IN
! 468: * MB = number of eigenvalues to compute
! 469: MB = IN
! 470: WEND = WBEGIN + MB - 1
! 471: * Define 1/4 and 3/4 points of the spectrum
! 472: S1 = ISLEFT + FOURTH * SPDIAM
! 473: S2 = ISRGHT - FOURTH * SPDIAM
! 474: ELSE
! 475: * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
! 476: * approximation.
! 477: * choose sigma
! 478: IF( USEDQD ) THEN
! 479: S1 = ISLEFT + FOURTH * SPDIAM
! 480: S2 = ISRGHT - FOURTH * SPDIAM
! 481: ELSE
! 482: TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL)
! 483: S1 = MAX(ISLEFT,VL) + FOURTH * TMP
! 484: S2 = MIN(ISRGHT,VU) - FOURTH * TMP
! 485: ENDIF
! 486: ENDIF
! 487:
! 488: * Compute the negcount at the 1/4 and 3/4 points
! 489: IF(MB.GT.1) THEN
! 490: CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN),
! 491: $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
! 492: ENDIF
! 493:
! 494: IF(MB.EQ.1) THEN
! 495: SIGMA = GL
! 496: SGNDEF = ONE
! 497: ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
! 498: IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
! 499: SIGMA = MAX(ISLEFT,GL)
! 500: ELSEIF( USEDQD ) THEN
! 501: * use Gerschgorin bound as shift to get pos def matrix
! 502: * for dqds
! 503: SIGMA = ISLEFT
! 504: ELSE
! 505: * use approximation of the first desired eigenvalue of the
! 506: * block as shift
! 507: SIGMA = MAX(ISLEFT,VL)
! 508: ENDIF
! 509: SGNDEF = ONE
! 510: ELSE
! 511: IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
! 512: SIGMA = MIN(ISRGHT,GU)
! 513: ELSEIF( USEDQD ) THEN
! 514: * use Gerschgorin bound as shift to get neg def matrix
! 515: * for dqds
! 516: SIGMA = ISRGHT
! 517: ELSE
! 518: * use approximation of the first desired eigenvalue of the
! 519: * block as shift
! 520: SIGMA = MIN(ISRGHT,VU)
! 521: ENDIF
! 522: SGNDEF = -ONE
! 523: ENDIF
! 524:
! 525:
! 526: * An initial SIGMA has been chosen that will be used for computing
! 527: * T - SIGMA I = L D L^T
! 528: * Define the increment TAU of the shift in case the initial shift
! 529: * needs to be refined to obtain a factorization with not too much
! 530: * element growth.
! 531: IF( USEDQD ) THEN
! 532: * The initial SIGMA was to the outer end of the spectrum
! 533: * the matrix is definite and we need not retreat.
! 534: TAU = SPDIAM*EPS*N + TWO*PIVMIN
! 535: ELSE
! 536: IF(MB.GT.1) THEN
! 537: CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
! 538: AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN))
! 539: IF( SGNDEF.EQ.ONE ) THEN
! 540: TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
! 541: TAU = MAX(TAU,WERR(WBEGIN))
! 542: ELSE
! 543: TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
! 544: TAU = MAX(TAU,WERR(WEND))
! 545: ENDIF
! 546: ELSE
! 547: TAU = WERR(WBEGIN)
! 548: ENDIF
! 549: ENDIF
! 550: *
! 551: DO 80 IDUM = 1, MAXTRY
! 552: * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
! 553: * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
! 554: * pivots in WORK(2*IN+1:3*IN)
! 555: DPIVOT = D( IBEGIN ) - SIGMA
! 556: WORK( 1 ) = DPIVOT
! 557: DMAX = ABS( WORK(1) )
! 558: J = IBEGIN
! 559: DO 70 I = 1, IN - 1
! 560: WORK( 2*IN+I ) = ONE / WORK( I )
! 561: TMP = E( J )*WORK( 2*IN+I )
! 562: WORK( IN+I ) = TMP
! 563: DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
! 564: WORK( I+1 ) = DPIVOT
! 565: DMAX = MAX( DMAX, ABS(DPIVOT) )
! 566: J = J + 1
! 567: 70 CONTINUE
! 568: * check for element growth
! 569: IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
! 570: NOREP = .TRUE.
! 571: ELSE
! 572: NOREP = .FALSE.
! 573: ENDIF
! 574: IF( USEDQD .AND. .NOT.NOREP ) THEN
! 575: * Ensure the definiteness of the representation
! 576: * All entries of D (of L D L^T) must have the same sign
! 577: DO 71 I = 1, IN
! 578: TMP = SGNDEF*WORK( I )
! 579: IF( TMP.LT.ZERO ) NOREP = .TRUE.
! 580: 71 CONTINUE
! 581: ENDIF
! 582: IF(NOREP) THEN
! 583: * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
! 584: * shift which makes the matrix definite. So we should end up
! 585: * here really only in the case of IRANGE = VALRNG or INDRNG.
! 586: IF( IDUM.EQ.MAXTRY-1 ) THEN
! 587: IF( SGNDEF.EQ.ONE ) THEN
! 588: * The fudged Gerschgorin shift should succeed
! 589: SIGMA =
! 590: $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
! 591: ELSE
! 592: SIGMA =
! 593: $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
! 594: END IF
! 595: ELSE
! 596: SIGMA = SIGMA - SGNDEF * TAU
! 597: TAU = TWO * TAU
! 598: END IF
! 599: ELSE
! 600: * an initial RRR is found
! 601: GO TO 83
! 602: END IF
! 603: 80 CONTINUE
! 604: * if the program reaches this point, no base representation could be
! 605: * found in MAXTRY iterations.
! 606: INFO = 2
! 607: RETURN
! 608:
! 609: 83 CONTINUE
! 610: * At this point, we have found an initial base representation
! 611: * T - SIGMA I = L D L^T with not too much element growth.
! 612: * Store the shift.
! 613: E( IEND ) = SIGMA
! 614: * Store D and L.
! 615: CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
! 616: CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
! 617:
! 618:
! 619: IF(MB.GT.1 ) THEN
! 620: *
! 621: * Perturb each entry of the base representation by a small
! 622: * (but random) relative amount to overcome difficulties with
! 623: * glued matrices.
! 624: *
! 625: DO 122 I = 1, 4
! 626: ISEED( I ) = 1
! 627: 122 CONTINUE
! 628:
! 629: CALL DLARNV(2, ISEED, 2*IN-1, WORK(1))
! 630: DO 125 I = 1,IN-1
! 631: D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
! 632: E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
! 633: 125 CONTINUE
! 634: D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
! 635: *
! 636: ENDIF
! 637: *
! 638: * Don't update the Gerschgorin intervals because keeping track
! 639: * of the updates would be too much work in DLARRV.
! 640: * We update W instead and use it to locate the proper Gerschgorin
! 641: * intervals.
! 642:
! 643: * Compute the required eigenvalues of L D L' by bisection or dqds
! 644: IF ( .NOT.USEDQD ) THEN
! 645: * If DLARRD has been used, shift the eigenvalue approximations
! 646: * according to their representation. This is necessary for
! 647: * a uniform DLARRV since dqds computes eigenvalues of the
! 648: * shifted representation. In DLARRV, W will always hold the
! 649: * UNshifted eigenvalue approximation.
! 650: DO 134 J=WBEGIN,WEND
! 651: W(J) = W(J) - SIGMA
! 652: WERR(J) = WERR(J) + ABS(W(J)) * EPS
! 653: 134 CONTINUE
! 654: * call DLARRB to reduce eigenvalue error of the approximations
! 655: * from DLARRD
! 656: DO 135 I = IBEGIN, IEND-1
! 657: WORK( I ) = D( I ) * E( I )**2
! 658: 135 CONTINUE
! 659: * use bisection to find EV from INDL to INDU
! 660: CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN),
! 661: $ INDL, INDU, RTOL1, RTOL2, INDL-1,
! 662: $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
! 663: $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
! 664: $ IN, IINFO )
! 665: IF( IINFO .NE. 0 ) THEN
! 666: INFO = -4
! 667: RETURN
! 668: END IF
! 669: * DLARRB computes all gaps correctly except for the last one
! 670: * Record distance to VU/GU
! 671: WGAP( WEND ) = MAX( ZERO,
! 672: $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
! 673: DO 138 I = INDL, INDU
! 674: M = M + 1
! 675: IBLOCK(M) = JBLK
! 676: INDEXW(M) = I
! 677: 138 CONTINUE
! 678: ELSE
! 679: * Call dqds to get all eigs (and then possibly delete unwanted
! 680: * eigenvalues).
! 681: * Note that dqds finds the eigenvalues of the L D L^T representation
! 682: * of T to high relative accuracy. High relative accuracy
! 683: * might be lost when the shift of the RRR is subtracted to obtain
! 684: * the eigenvalues of T. However, T is not guaranteed to define its
! 685: * eigenvalues to high relative accuracy anyway.
! 686: * Set RTOL to the order of the tolerance used in DLASQ2
! 687: * This is an ESTIMATED error, the worst case bound is 4*N*EPS
! 688: * which is usually too large and requires unnecessary work to be
! 689: * done by bisection when computing the eigenvectors
! 690: RTOL = LOG(DBLE(IN)) * FOUR * EPS
! 691: J = IBEGIN
! 692: DO 140 I = 1, IN - 1
! 693: WORK( 2*I-1 ) = ABS( D( J ) )
! 694: WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
! 695: J = J + 1
! 696: 140 CONTINUE
! 697: WORK( 2*IN-1 ) = ABS( D( IEND ) )
! 698: WORK( 2*IN ) = ZERO
! 699: CALL DLASQ2( IN, WORK, IINFO )
! 700: IF( IINFO .NE. 0 ) THEN
! 701: * If IINFO = -5 then an index is part of a tight cluster
! 702: * and should be changed. The index is in IWORK(1) and the
! 703: * gap is in WORK(N+1)
! 704: INFO = -5
! 705: RETURN
! 706: ELSE
! 707: * Test that all eigenvalues are positive as expected
! 708: DO 149 I = 1, IN
! 709: IF( WORK( I ).LT.ZERO ) THEN
! 710: INFO = -6
! 711: RETURN
! 712: ENDIF
! 713: 149 CONTINUE
! 714: END IF
! 715: IF( SGNDEF.GT.ZERO ) THEN
! 716: DO 150 I = INDL, INDU
! 717: M = M + 1
! 718: W( M ) = WORK( IN-I+1 )
! 719: IBLOCK( M ) = JBLK
! 720: INDEXW( M ) = I
! 721: 150 CONTINUE
! 722: ELSE
! 723: DO 160 I = INDL, INDU
! 724: M = M + 1
! 725: W( M ) = -WORK( I )
! 726: IBLOCK( M ) = JBLK
! 727: INDEXW( M ) = I
! 728: 160 CONTINUE
! 729: END IF
! 730:
! 731: DO 165 I = M - MB + 1, M
! 732: * the value of RTOL below should be the tolerance in DLASQ2
! 733: WERR( I ) = RTOL * ABS( W(I) )
! 734: 165 CONTINUE
! 735: DO 166 I = M - MB + 1, M - 1
! 736: * compute the right gap between the intervals
! 737: WGAP( I ) = MAX( ZERO,
! 738: $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
! 739: 166 CONTINUE
! 740: WGAP( M ) = MAX( ZERO,
! 741: $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
! 742: END IF
! 743: * proceed with next block
! 744: IBEGIN = IEND + 1
! 745: WBEGIN = WEND + 1
! 746: 170 CONTINUE
! 747: *
! 748:
! 749: RETURN
! 750: *
! 751: * end of DLARRE
! 752: *
! 753: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarre.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack