Annotation of rpl/lapack/lapack/dlarrk.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLARRK( N, IW, GL, GU,
! 2: $ D, E2, PIVMIN, RELTOL, W, WERR, INFO)
! 3: IMPLICIT NONE
! 4: *
! 5: * -- LAPACK auxiliary 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: *
! 10: * .. Scalar Arguments ..
! 11: INTEGER INFO, IW, N
! 12: DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR
! 13: * ..
! 14: * .. Array Arguments ..
! 15: DOUBLE PRECISION D( * ), E2( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DLARRK computes one eigenvalue of a symmetric tridiagonal
! 22: * matrix T to suitable accuracy. This is an auxiliary code to be
! 23: * called from DSTEMR.
! 24: *
! 25: * To avoid overflow, the matrix must be scaled so that its
! 26: * largest element is no greater than overflow**(1/2) *
! 27: * underflow**(1/4) in absolute value, and for greatest
! 28: * accuracy, it should not be much smaller than that.
! 29: *
! 30: * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
! 31: * Matrix", Report CS41, Computer Science Dept., Stanford
! 32: * University, July 21, 1966.
! 33: *
! 34: * Arguments
! 35: * =========
! 36: *
! 37: * N (input) INTEGER
! 38: * The order of the tridiagonal matrix T. N >= 0.
! 39: *
! 40: * IW (input) INTEGER
! 41: * The index of the eigenvalues to be returned.
! 42: *
! 43: * GL (input) DOUBLE PRECISION
! 44: * GU (input) DOUBLE PRECISION
! 45: * An upper and a lower bound on the eigenvalue.
! 46: *
! 47: * D (input) DOUBLE PRECISION array, dimension (N)
! 48: * The n diagonal elements of the tridiagonal matrix T.
! 49: *
! 50: * E2 (input) DOUBLE PRECISION array, dimension (N-1)
! 51: * The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
! 52: *
! 53: * PIVMIN (input) DOUBLE PRECISION
! 54: * The minimum pivot allowed in the Sturm sequence for T.
! 55: *
! 56: * RELTOL (input) DOUBLE PRECISION
! 57: * The minimum relative width of an interval. When an interval
! 58: * is narrower than RELTOL times the larger (in
! 59: * magnitude) endpoint, then it is considered to be
! 60: * sufficiently small, i.e., converged. Note: this should
! 61: * always be at least radix*machine epsilon.
! 62: *
! 63: * W (output) DOUBLE PRECISION
! 64: *
! 65: * WERR (output) DOUBLE PRECISION
! 66: * The error bound on the corresponding eigenvalue approximation
! 67: * in W.
! 68: *
! 69: * INFO (output) INTEGER
! 70: * = 0: Eigenvalue converged
! 71: * = -1: Eigenvalue did NOT converge
! 72: *
! 73: * Internal Parameters
! 74: * ===================
! 75: *
! 76: * FUDGE DOUBLE PRECISION, default = 2
! 77: * A "fudge factor" to widen the Gershgorin intervals.
! 78: *
! 79: * =====================================================================
! 80: *
! 81: * .. Parameters ..
! 82: DOUBLE PRECISION FUDGE, HALF, TWO, ZERO
! 83: PARAMETER ( HALF = 0.5D0, TWO = 2.0D0,
! 84: $ FUDGE = TWO, ZERO = 0.0D0 )
! 85: * ..
! 86: * .. Local Scalars ..
! 87: INTEGER I, IT, ITMAX, NEGCNT
! 88: DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
! 89: $ TMP2, TNORM
! 90: * ..
! 91: * .. External Functions ..
! 92: DOUBLE PRECISION DLAMCH
! 93: EXTERNAL DLAMCH
! 94: * ..
! 95: * .. Intrinsic Functions ..
! 96: INTRINSIC ABS, INT, LOG, MAX
! 97: * ..
! 98: * .. Executable Statements ..
! 99: *
! 100: * Get machine constants
! 101: EPS = DLAMCH( 'P' )
! 102:
! 103: TNORM = MAX( ABS( GL ), ABS( GU ) )
! 104: RTOLI = RELTOL
! 105: ATOLI = FUDGE*TWO*PIVMIN
! 106:
! 107: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
! 108: $ LOG( TWO ) ) + 2
! 109:
! 110: INFO = -1
! 111:
! 112: LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
! 113: RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
! 114: IT = 0
! 115:
! 116: 10 CONTINUE
! 117: *
! 118: * Check if interval converged or maximum number of iterations reached
! 119: *
! 120: TMP1 = ABS( RIGHT - LEFT )
! 121: TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
! 122: IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
! 123: INFO = 0
! 124: GOTO 30
! 125: ENDIF
! 126: IF(IT.GT.ITMAX)
! 127: $ GOTO 30
! 128:
! 129: *
! 130: * Count number of negative pivots for mid-point
! 131: *
! 132: IT = IT + 1
! 133: MID = HALF * (LEFT + RIGHT)
! 134: NEGCNT = 0
! 135: TMP1 = D( 1 ) - MID
! 136: IF( ABS( TMP1 ).LT.PIVMIN )
! 137: $ TMP1 = -PIVMIN
! 138: IF( TMP1.LE.ZERO )
! 139: $ NEGCNT = NEGCNT + 1
! 140: *
! 141: DO 20 I = 2, N
! 142: TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
! 143: IF( ABS( TMP1 ).LT.PIVMIN )
! 144: $ TMP1 = -PIVMIN
! 145: IF( TMP1.LE.ZERO )
! 146: $ NEGCNT = NEGCNT + 1
! 147: 20 CONTINUE
! 148:
! 149: IF(NEGCNT.GE.IW) THEN
! 150: RIGHT = MID
! 151: ELSE
! 152: LEFT = MID
! 153: ENDIF
! 154: GOTO 10
! 155:
! 156: 30 CONTINUE
! 157: *
! 158: * Converged or maximum number of iterations reached
! 159: *
! 160: W = HALF * (LEFT + RIGHT)
! 161: WERR = HALF * ABS( RIGHT - LEFT )
! 162:
! 163: RETURN
! 164: *
! 165: * End of DLARRK
! 166: *
! 167: END
CVSweb interface <joel.bertrand@systella.fr>