1*> \brief \b DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLARRK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLARRK( N, IW, GL, GU,
22*                           D, E2, PIVMIN, RELTOL, W, WERR, INFO)
23*
24*       .. Scalar Arguments ..
25*       INTEGER   INFO, IW, N
26*       DOUBLE PRECISION    PIVMIN, RELTOL, GL, GU, W, WERR
27*       ..
28*       .. Array Arguments ..
29*       DOUBLE PRECISION   D( * ), E2( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> DLARRK computes one eigenvalue of a symmetric tridiagonal
39*> matrix T to suitable accuracy. This is an auxiliary code to be
40*> called from DSTEMR.
41*>
42*> To avoid overflow, the matrix must be scaled so that its
43*> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
44*> accuracy, it should not be much smaller than that.
45*>
46*> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
47*> Matrix", Report CS41, Computer Science Dept., Stanford
48*> University, July 21, 1966.
49*> \endverbatim
50*
51*  Arguments:
52*  ==========
53*
54*> \param[in] N
55*> \verbatim
56*>          N is INTEGER
57*>          The order of the tridiagonal matrix T.  N >= 0.
58*> \endverbatim
59*>
60*> \param[in] IW
61*> \verbatim
62*>          IW is INTEGER
63*>          The index of the eigenvalues to be returned.
64*> \endverbatim
65*>
66*> \param[in] GL
67*> \verbatim
68*>          GL is DOUBLE PRECISION
69*> \endverbatim
70*>
71*> \param[in] GU
72*> \verbatim
73*>          GU is DOUBLE PRECISION
74*>          An upper and a lower bound on the eigenvalue.
75*> \endverbatim
76*>
77*> \param[in] D
78*> \verbatim
79*>          D is DOUBLE PRECISION array, dimension (N)
80*>          The n diagonal elements of the tridiagonal matrix T.
81*> \endverbatim
82*>
83*> \param[in] E2
84*> \verbatim
85*>          E2 is DOUBLE PRECISION array, dimension (N-1)
86*>          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
87*> \endverbatim
88*>
89*> \param[in] PIVMIN
90*> \verbatim
91*>          PIVMIN is DOUBLE PRECISION
92*>          The minimum pivot allowed in the Sturm sequence for T.
93*> \endverbatim
94*>
95*> \param[in] RELTOL
96*> \verbatim
97*>          RELTOL is DOUBLE PRECISION
98*>          The minimum relative width of an interval.  When an interval
99*>          is narrower than RELTOL times the larger (in
100*>          magnitude) endpoint, then it is considered to be
101*>          sufficiently small, i.e., converged.  Note: this should
102*>          always be at least radix*machine epsilon.
103*> \endverbatim
104*>
105*> \param[out] W
106*> \verbatim
107*>          W is DOUBLE PRECISION
108*> \endverbatim
109*>
110*> \param[out] WERR
111*> \verbatim
112*>          WERR is DOUBLE PRECISION
113*>          The error bound on the corresponding eigenvalue approximation
114*>          in W.
115*> \endverbatim
116*>
117*> \param[out] INFO
118*> \verbatim
119*>          INFO is INTEGER
120*>          = 0:       Eigenvalue converged
121*>          = -1:      Eigenvalue did NOT converge
122*> \endverbatim
123*
124*> \par Internal Parameters:
125*  =========================
126*>
127*> \verbatim
128*>  FUDGE   DOUBLE PRECISION, default = 2
129*>          A "fudge factor" to widen the Gershgorin intervals.
130*> \endverbatim
131*
132*  Authors:
133*  ========
134*
135*> \author Univ. of Tennessee
136*> \author Univ. of California Berkeley
137*> \author Univ. of Colorado Denver
138*> \author NAG Ltd.
139*
140*> \ingroup OTHERauxiliary
141*
142*  =====================================================================
143      SUBROUTINE DLARRK( N, IW, GL, GU,
144     $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
145*
146*  -- LAPACK auxiliary routine --
147*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
148*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149*
150*     .. Scalar Arguments ..
151      INTEGER   INFO, IW, N
152      DOUBLE PRECISION    PIVMIN, RELTOL, GL, GU, W, WERR
153*     ..
154*     .. Array Arguments ..
155      DOUBLE PRECISION   D( * ), E2( * )
156*     ..
157*
158*  =====================================================================
159*
160*     .. Parameters ..
161      DOUBLE PRECISION   FUDGE, HALF, TWO, ZERO
162      PARAMETER          ( HALF = 0.5D0, TWO = 2.0D0,
163     $                     FUDGE = TWO, ZERO = 0.0D0 )
164*     ..
165*     .. Local Scalars ..
166      INTEGER   I, IT, ITMAX, NEGCNT
167      DOUBLE PRECISION   ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
168     $                   TMP2, TNORM
169*     ..
170*     .. External Functions ..
171      DOUBLE PRECISION   DLAMCH
172      EXTERNAL   DLAMCH
173*     ..
174*     .. Intrinsic Functions ..
175      INTRINSIC          ABS, INT, LOG, MAX
176*     ..
177*     .. Executable Statements ..
178*
179*     Quick return if possible
180*
181      IF( N.LE.0 ) THEN
182         INFO = 0
183         RETURN
184      END IF
185*
186*     Get machine constants
187      EPS = DLAMCH( 'P' )
188
189      TNORM = MAX( ABS( GL ), ABS( GU ) )
190      RTOLI = RELTOL
191      ATOLI = FUDGE*TWO*PIVMIN
192
193      ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
194     $           LOG( TWO ) ) + 2
195
196      INFO = -1
197
198      LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
199      RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
200      IT = 0
201
202 10   CONTINUE
203*
204*     Check if interval converged or maximum number of iterations reached
205*
206      TMP1 = ABS( RIGHT - LEFT )
207      TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
208      IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
209         INFO = 0
210         GOTO 30
211      ENDIF
212      IF(IT.GT.ITMAX)
213     $   GOTO 30
214
215*
216*     Count number of negative pivots for mid-point
217*
218      IT = IT + 1
219      MID = HALF * (LEFT + RIGHT)
220      NEGCNT = 0
221      TMP1 = D( 1 ) - MID
222      IF( ABS( TMP1 ).LT.PIVMIN )
223     $   TMP1 = -PIVMIN
224      IF( TMP1.LE.ZERO )
225     $   NEGCNT = NEGCNT + 1
226*
227      DO 20 I = 2, N
228         TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
229         IF( ABS( TMP1 ).LT.PIVMIN )
230     $      TMP1 = -PIVMIN
231         IF( TMP1.LE.ZERO )
232     $      NEGCNT = NEGCNT + 1
233 20   CONTINUE
234
235      IF(NEGCNT.GE.IW) THEN
236         RIGHT = MID
237      ELSE
238         LEFT = MID
239      ENDIF
240      GOTO 10
241
242 30   CONTINUE
243*
244*     Converged or maximum number of iterations reached
245*
246      W = HALF * (LEFT + RIGHT)
247      WERR = HALF * ABS( RIGHT - LEFT )
248
249      RETURN
250*
251*     End of DLARRK
252*
253      END
254