001:       SUBROUTINE SLARRK( N, IW, GL, GU,
002:      $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
003:       IMPLICIT NONE
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       INTEGER   INFO, IW, N
012:       REAL                PIVMIN, RELTOL, GL, GU, W, WERR
013: *     ..
014: *     .. Array Arguments ..
015:       REAL               D( * ), E2( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SLARRK computes one eigenvalue of a symmetric tridiagonal
022: *  matrix T to suitable accuracy. This is an auxiliary code to be
023: *  called from SSTEMR.
024: *
025: *  To avoid overflow, the matrix must be scaled so that its
026: *  largest element is no greater than overflow**(1/2) *
027: *  underflow**(1/4) in absolute value, and for greatest
028: *  accuracy, it should not be much smaller than that.
029: *
030: *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
031: *  Matrix", Report CS41, Computer Science Dept., Stanford
032: *  University, July 21, 1966.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  N       (input) INTEGER
038: *          The order of the tridiagonal matrix T.  N >= 0.
039: *
040: *  IW      (input) INTEGER
041: *          The index of the eigenvalues to be returned.
042: *
043: *  GL      (input) REAL            
044: *  GU      (input) REAL            
045: *          An upper and a lower bound on the eigenvalue.
046: *
047: *  D       (input) REAL             array, dimension (N)
048: *          The n diagonal elements of the tridiagonal matrix T.
049: *
050: *  E2      (input) REAL             array, dimension (N-1)
051: *          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
052: *
053: *  PIVMIN  (input) REAL            
054: *          The minimum pivot allowed in the Sturm sequence for T.
055: *
056: *  RELTOL  (input) REAL            
057: *          The minimum relative width of an interval.  When an interval
058: *          is narrower than RELTOL times the larger (in
059: *          magnitude) endpoint, then it is considered to be
060: *          sufficiently small, i.e., converged.  Note: this should
061: *          always be at least radix*machine epsilon.
062: *
063: *  W       (output) REAL            
064: *
065: *  WERR    (output) REAL            
066: *          The error bound on the corresponding eigenvalue approximation
067: *          in W.
068: *
069: *  INFO    (output) INTEGER
070: *          = 0:       Eigenvalue converged
071: *          = -1:      Eigenvalue did NOT converge
072: *
073: *  Internal Parameters
074: *  ===================
075: *
076: *  FUDGE   REAL            , default = 2
077: *          A "fudge factor" to widen the Gershgorin intervals.
078: *
079: *  =====================================================================
080: *
081: *     .. Parameters ..
082:       REAL               FUDGE, HALF, TWO, ZERO
083:       PARAMETER          ( HALF = 0.5E0, TWO = 2.0E0,
084:      $                     FUDGE = TWO, ZERO = 0.0E0 )
085: *     ..
086: *     .. Local Scalars ..
087:       INTEGER   I, IT, ITMAX, NEGCNT
088:       REAL               ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
089:      $                   TMP2, TNORM
090: *     ..
091: *     .. External Functions ..
092:       REAL               SLAMCH
093:       EXTERNAL   SLAMCH
094: *     ..
095: *     .. Intrinsic Functions ..
096:       INTRINSIC          ABS, INT, LOG, MAX
097: *     ..
098: *     .. Executable Statements ..
099: *
100: *     Get machine constants
101:       EPS = SLAMCH( '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 SLARRK
166: *
167:       END
168: