LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \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