LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlarrk.f
Go to the documentation of this file.
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
subroutine dlarrk(N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
Definition: dlarrk.f:145