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*> \date September 2012
141*
142*> \ingroup auxOTHERauxiliary
143*
144*  =====================================================================
145      SUBROUTINE DLARRK( N, IW, GL, GU,
146     $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
147*
148*  -- LAPACK auxiliary routine (version 3.4.2) --
149*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
150*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151*     September 2012
152*
153*     .. Scalar Arguments ..
154      INTEGER   INFO, IW, N
155      DOUBLE PRECISION    PIVMIN, RELTOL, GL, GU, W, WERR
156*     ..
157*     .. Array Arguments ..
158      DOUBLE PRECISION   D( * ), E2( * )
159*     ..
160*
161*  =====================================================================
162*
163*     .. Parameters ..
164      DOUBLE PRECISION   FUDGE, HALF, TWO, ZERO
165      PARAMETER          ( HALF = 0.5D0, TWO = 2.0D0,
166     $                     FUDGE = TWO, ZERO = 0.0D0 )
167*     ..
168*     .. Local Scalars ..
169      INTEGER   I, IT, ITMAX, NEGCNT
170      DOUBLE PRECISION   ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
171     $                   TMP2, TNORM
172*     ..
173*     .. External Functions ..
174      DOUBLE PRECISION   DLAMCH
175      EXTERNAL   DLAMCH
176*     ..
177*     .. Intrinsic Functions ..
178      INTRINSIC          ABS, INT, LOG, MAX
179*     ..
180*     .. Executable Statements ..
181*
182*     Get machine constants
183      EPS = DLAMCH( 'P' )
184
185      TNORM = MAX( ABS( GL ), ABS( GU ) )
186      RTOLI = RELTOL
187      ATOLI = FUDGE*TWO*PIVMIN
188
189      ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
190     $           LOG( TWO ) ) + 2
191
192      INFO = -1
193
194      LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
195      RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
196      IT = 0
197
198 10   CONTINUE
199*
200*     Check if interval converged or maximum number of iterations reached
201*
202      TMP1 = ABS( RIGHT - LEFT )
203      TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
204      IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
205         INFO = 0
206         GOTO 30
207      ENDIF
208      IF(IT.GT.ITMAX)
209     $   GOTO 30
210
211*
212*     Count number of negative pivots for mid-point
213*
214      IT = IT + 1
215      MID = HALF * (LEFT + RIGHT)
216      NEGCNT = 0
217      TMP1 = D( 1 ) - MID
218      IF( ABS( TMP1 ).LT.PIVMIN )
219     $   TMP1 = -PIVMIN
220      IF( TMP1.LE.ZERO )
221     $   NEGCNT = NEGCNT + 1
222*
223      DO 20 I = 2, N
224         TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
225         IF( ABS( TMP1 ).LT.PIVMIN )
226     $      TMP1 = -PIVMIN
227         IF( TMP1.LE.ZERO )
228     $      NEGCNT = NEGCNT + 1
229 20   CONTINUE
230
231      IF(NEGCNT.GE.IW) THEN
232         RIGHT = MID
233      ELSE
234         LEFT = MID
235      ENDIF
236      GOTO 10
237
238 30   CONTINUE
239*
240*     Converged or maximum number of iterations reached
241*
242      W = HALF * (LEFT + RIGHT)
243      WERR = HALF * ABS( RIGHT - LEFT )
244
245      RETURN
246*
247*     End of DLARRK
248*
249      END
250