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