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: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER   INFO, IW, N
011:       REAL                PIVMIN, RELTOL, GL, GU, W, WERR
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               D( * ), E2( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SLARRK computes one eigenvalue of a symmetric tridiagonal
021: *  matrix T to suitable accuracy. This is an auxiliary code to be
022: *  called from SSTEMR.
023: *
024: *  To avoid overflow, the matrix must be scaled so that its
025: *  largest element is no greater than overflow**(1/2) *
026: *  underflow**(1/4) in absolute value, and for greatest
027: *  accuracy, it should not be much smaller than that.
028: *
029: *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
030: *  Matrix", Report CS41, Computer Science Dept., Stanford
031: *  University, July 21, 1966.
032: *
033: *  Arguments
034: *  =========
035: *
036: *  N       (input) INTEGER
037: *          The order of the tridiagonal matrix T.  N >= 0.
038: *
039: *  IW      (input) INTEGER
040: *          The index of the eigenvalues to be returned.
041: *
042: *  GL      (input) REAL            
043: *  GU      (input) REAL            
044: *          An upper and a lower bound on the eigenvalue.
045: *
046: *  D       (input) REAL             array, dimension (N)
047: *          The n diagonal elements of the tridiagonal matrix T.
048: *
049: *  E2      (input) REAL             array, dimension (N-1)
050: *          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
051: *
052: *  PIVMIN  (input) REAL            
053: *          The minimum pivot allowed in the Sturm sequence for T.
054: *
055: *  RELTOL  (input) REAL            
056: *          The minimum relative width of an interval.  When an interval
057: *          is narrower than RELTOL times the larger (in
058: *          magnitude) endpoint, then it is considered to be
059: *          sufficiently small, i.e., converged.  Note: this should
060: *          always be at least radix*machine epsilon.
061: *
062: *  W       (output) REAL            
063: *
064: *  WERR    (output) REAL            
065: *          The error bound on the corresponding eigenvalue approximation
066: *          in W.
067: *
068: *  INFO    (output) INTEGER
069: *          = 0:       Eigenvalue converged
070: *          = -1:      Eigenvalue did NOT converge
071: *
072: *  Internal Parameters
073: *  ===================
074: *
075: *  FUDGE   REAL            , default = 2
076: *          A "fudge factor" to widen the Gershgorin intervals.
077: *
078: *  =====================================================================
079: *
080: *     .. Parameters ..
081:       REAL               FUDGE, HALF, TWO, ZERO
082:       PARAMETER          ( HALF = 0.5E0, TWO = 2.0E0,
083:      $                     FUDGE = TWO, ZERO = 0.0E0 )
084: *     ..
085: *     .. Local Scalars ..
086:       INTEGER   I, IT, ITMAX, NEGCNT
087:       REAL               ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
088:      $                   TMP2, TNORM
089: *     ..
090: *     .. External Functions ..
091:       REAL               SLAMCH
092:       EXTERNAL   SLAMCH
093: *     ..
094: *     .. Intrinsic Functions ..
095:       INTRINSIC          ABS, INT, LOG, MAX
096: *     ..
097: *     .. Executable Statements ..
098: *
099: *     Get machine constants
100:       EPS = SLAMCH( 'P' )
101: 
102:       TNORM = MAX( ABS( GL ), ABS( GU ) )
103:       RTOLI = RELTOL
104:       ATOLI = FUDGE*TWO*PIVMIN
105: 
106:       ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
107:      $           LOG( TWO ) ) + 2
108: 
109:       INFO = -1
110: 
111:       LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
112:       RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
113:       IT = 0
114: 
115:  10   CONTINUE
116: *
117: *     Check if interval converged or maximum number of iterations reached
118: *
119:       TMP1 = ABS( RIGHT - LEFT )
120:       TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
121:       IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
122:          INFO = 0
123:          GOTO 30
124:       ENDIF
125:       IF(IT.GT.ITMAX)
126:      $   GOTO 30
127: 
128: *
129: *     Count number of negative pivots for mid-point
130: *
131:       IT = IT + 1
132:       MID = HALF * (LEFT + RIGHT)
133:       NEGCNT = 0
134:       TMP1 = D( 1 ) - MID
135:       IF( ABS( TMP1 ).LT.PIVMIN )
136:      $   TMP1 = -PIVMIN
137:       IF( TMP1.LE.ZERO )
138:      $   NEGCNT = NEGCNT + 1
139: *
140:       DO 20 I = 2, N
141:          TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
142:          IF( ABS( TMP1 ).LT.PIVMIN )
143:      $      TMP1 = -PIVMIN
144:          IF( TMP1.LE.ZERO )
145:      $      NEGCNT = NEGCNT + 1
146:  20   CONTINUE
147: 
148:       IF(NEGCNT.GE.IW) THEN
149:          RIGHT = MID
150:       ELSE
151:          LEFT = MID
152:       ENDIF
153:       GOTO 10
154: 
155:  30   CONTINUE
156: *
157: *     Converged or maximum number of iterations reached
158: *
159:       W = HALF * (LEFT + RIGHT)
160:       WERR = HALF * ABS( RIGHT - LEFT )
161: 
162:       RETURN
163: *
164: *     End of SLARRK
165: *
166:       END
167: