LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlarrj.f
Go to the documentation of this file.
1*> \brief \b DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DLARRJ + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrj.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrj.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrj.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
20* RTOL, OFFSET, W, WERR, WORK, IWORK,
21* PIVMIN, SPDIAM, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER IFIRST, ILAST, INFO, N, OFFSET
25* DOUBLE PRECISION PIVMIN, RTOL, SPDIAM
26* ..
27* .. Array Arguments ..
28* INTEGER IWORK( * )
29* DOUBLE PRECISION D( * ), E2( * ), W( * ),
30* $ WERR( * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> Given the initial eigenvalue approximations of T, DLARRJ
40*> does bisection to refine the eigenvalues of T,
41*> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
42*> guesses for these eigenvalues are input in W, the corresponding estimate
43*> of the error in these guesses in WERR. During bisection, intervals
44*> [left, right] are maintained by storing their mid-points and
45*> semi-widths in the arrays W and WERR respectively.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] N
52*> \verbatim
53*> N is INTEGER
54*> The order of the matrix.
55*> \endverbatim
56*>
57*> \param[in] D
58*> \verbatim
59*> D is DOUBLE PRECISION array, dimension (N)
60*> The N diagonal elements of T.
61*> \endverbatim
62*>
63*> \param[in] E2
64*> \verbatim
65*> E2 is DOUBLE PRECISION array, dimension (N-1)
66*> The Squares of the (N-1) subdiagonal elements of T.
67*> \endverbatim
68*>
69*> \param[in] IFIRST
70*> \verbatim
71*> IFIRST is INTEGER
72*> The index of the first eigenvalue to be computed.
73*> \endverbatim
74*>
75*> \param[in] ILAST
76*> \verbatim
77*> ILAST is INTEGER
78*> The index of the last eigenvalue to be computed.
79*> \endverbatim
80*>
81*> \param[in] RTOL
82*> \verbatim
83*> RTOL is DOUBLE PRECISION
84*> Tolerance for the convergence of the bisection intervals.
85*> An interval [LEFT,RIGHT] has converged if
86*> RIGHT-LEFT < RTOL*MAX(|LEFT|,|RIGHT|).
87*> \endverbatim
88*>
89*> \param[in] OFFSET
90*> \verbatim
91*> OFFSET is INTEGER
92*> Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
93*> through ILAST-OFFSET elements of these arrays are to be used.
94*> \endverbatim
95*>
96*> \param[in,out] W
97*> \verbatim
98*> W is DOUBLE PRECISION array, dimension (N)
99*> On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
100*> estimates of the eigenvalues of L D L^T indexed IFIRST through
101*> ILAST.
102*> On output, these estimates are refined.
103*> \endverbatim
104*>
105*> \param[in,out] WERR
106*> \verbatim
107*> WERR is DOUBLE PRECISION array, dimension (N)
108*> On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
109*> the errors in the estimates of the corresponding elements in W.
110*> On output, these errors are refined.
111*> \endverbatim
112*>
113*> \param[out] WORK
114*> \verbatim
115*> WORK is DOUBLE PRECISION array, dimension (2*N)
116*> Workspace.
117*> \endverbatim
118*>
119*> \param[out] IWORK
120*> \verbatim
121*> IWORK is INTEGER array, dimension (2*N)
122*> Workspace.
123*> \endverbatim
124*>
125*> \param[in] PIVMIN
126*> \verbatim
127*> PIVMIN is DOUBLE PRECISION
128*> The minimum pivot in the Sturm sequence for T.
129*> \endverbatim
130*>
131*> \param[in] SPDIAM
132*> \verbatim
133*> SPDIAM is DOUBLE PRECISION
134*> The spectral diameter of T.
135*> \endverbatim
136*>
137*> \param[out] INFO
138*> \verbatim
139*> INFO is INTEGER
140*> Error flag.
141*> \endverbatim
142*
143* Authors:
144* ========
145*
146*> \author Univ. of Tennessee
147*> \author Univ. of California Berkeley
148*> \author Univ. of Colorado Denver
149*> \author NAG Ltd.
150*
151*> \ingroup larrj
152*
153*> \par Contributors:
154* ==================
155*>
156*> Beresford Parlett, University of California, Berkeley, USA \n
157*> Jim Demmel, University of California, Berkeley, USA \n
158*> Inderjit Dhillon, University of Texas, Austin, USA \n
159*> Osni Marques, LBNL/NERSC, USA \n
160*> Christof Voemel, University of California, Berkeley, USA
161*
162* =====================================================================
163 SUBROUTINE dlarrj( N, D, E2, IFIRST, ILAST,
164 $ RTOL, OFFSET, W, WERR, WORK, IWORK,
165 $ PIVMIN, SPDIAM, INFO )
166*
167* -- LAPACK auxiliary routine --
168* -- LAPACK is a software package provided by Univ. of Tennessee, --
169* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170*
171* .. Scalar Arguments ..
172 INTEGER IFIRST, ILAST, INFO, N, OFFSET
173 DOUBLE PRECISION PIVMIN, RTOL, SPDIAM
174* ..
175* .. Array Arguments ..
176 INTEGER IWORK( * )
177 DOUBLE PRECISION D( * ), E2( * ), W( * ),
178 $ werr( * ), work( * )
179* ..
180*
181* =====================================================================
182*
183* .. Parameters ..
184 DOUBLE PRECISION ZERO, ONE, TWO, HALF
185 PARAMETER ( ZERO = 0.0d0, one = 1.0d0, two = 2.0d0,
186 $ half = 0.5d0 )
187 INTEGER MAXITR
188* ..
189* .. Local Scalars ..
190 INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
191 $ OLNINT, P, PREV, SAVI1
192 DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
193*
194* ..
195* .. Intrinsic Functions ..
196 INTRINSIC abs, max
197* ..
198* .. Executable Statements ..
199*
200 info = 0
201*
202* Quick return if possible
203*
204 IF( n.LE.0 ) THEN
205 RETURN
206 END IF
207*
208 maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
209 $ log( two ) ) + 2
210*
211* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
212* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
213* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
214* for an unconverged interval is set to the index of the next unconverged
215* interval, and is -1 or 0 for a converged interval. Thus a linked
216* list of unconverged intervals is set up.
217*
218
219 i1 = ifirst
220 i2 = ilast
221* The number of unconverged intervals
222 nint = 0
223* The last unconverged interval found
224 prev = 0
225 DO 75 i = i1, i2
226 k = 2*i
227 ii = i - offset
228 left = w( ii ) - werr( ii )
229 mid = w(ii)
230 right = w( ii ) + werr( ii )
231 width = right - mid
232 tmp = max( abs( left ), abs( right ) )
233
234* The following test prevents the test of converged intervals
235 IF( width.LT.rtol*tmp ) THEN
236* This interval has already converged and does not need refinement.
237* (Note that the gaps might change through refining the
238* eigenvalues, however, they can only get bigger.)
239* Remove it from the list.
240 iwork( k-1 ) = -1
241* Make sure that I1 always points to the first unconverged interval
242 IF((i.EQ.i1).AND.(i.LT.i2)) i1 = i + 1
243 IF((prev.GE.i1).AND.(i.LE.i2)) iwork( 2*prev-1 ) = i + 1
244 ELSE
245* unconverged interval found
246 prev = i
247* Make sure that [LEFT,RIGHT] contains the desired eigenvalue
248*
249* Do while( CNT(LEFT).GT.I-1 )
250*
251 fac = one
252 20 CONTINUE
253 cnt = 0
254 s = left
255 dplus = d( 1 ) - s
256 IF( dplus.LT.zero ) cnt = cnt + 1
257 DO 30 j = 2, n
258 dplus = d( j ) - s - e2( j-1 )/dplus
259 IF( dplus.LT.zero ) cnt = cnt + 1
260 30 CONTINUE
261 IF( cnt.GT.i-1 ) THEN
262 left = left - werr( ii )*fac
263 fac = two*fac
264 GO TO 20
265 END IF
266*
267* Do while( CNT(RIGHT).LT.I )
268*
269 fac = one
270 50 CONTINUE
271 cnt = 0
272 s = right
273 dplus = d( 1 ) - s
274 IF( dplus.LT.zero ) cnt = cnt + 1
275 DO 60 j = 2, n
276 dplus = d( j ) - s - e2( j-1 )/dplus
277 IF( dplus.LT.zero ) cnt = cnt + 1
278 60 CONTINUE
279 IF( cnt.LT.i ) THEN
280 right = right + werr( ii )*fac
281 fac = two*fac
282 GO TO 50
283 END IF
284 nint = nint + 1
285 iwork( k-1 ) = i + 1
286 iwork( k ) = cnt
287 END IF
288 work( k-1 ) = left
289 work( k ) = right
290 75 CONTINUE
291
292
293 savi1 = i1
294*
295* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
296* and while (ITER.LT.MAXITR)
297*
298 iter = 0
299 80 CONTINUE
300 prev = i1 - 1
301 i = i1
302 olnint = nint
303
304 DO 100 p = 1, olnint
305 k = 2*i
306 ii = i - offset
307 next = iwork( k-1 )
308 left = work( k-1 )
309 right = work( k )
310 mid = half*( left + right )
311
312* semiwidth of interval
313 width = right - mid
314 tmp = max( abs( left ), abs( right ) )
315
316 IF( ( width.LT.rtol*tmp ) .OR.
317 $ (iter.EQ.maxitr) )THEN
318* reduce number of unconverged intervals
319 nint = nint - 1
320* Mark interval as converged.
321 iwork( k-1 ) = 0
322 IF( i1.EQ.i ) THEN
323 i1 = next
324 ELSE
325* Prev holds the last unconverged interval previously examined
326 IF(prev.GE.i1) iwork( 2*prev-1 ) = next
327 END IF
328 i = next
329 GO TO 100
330 END IF
331 prev = i
332*
333* Perform one bisection step
334*
335 cnt = 0
336 s = mid
337 dplus = d( 1 ) - s
338 IF( dplus.LT.zero ) cnt = cnt + 1
339 DO 90 j = 2, n
340 dplus = d( j ) - s - e2( j-1 )/dplus
341 IF( dplus.LT.zero ) cnt = cnt + 1
342 90 CONTINUE
343 IF( cnt.LE.i-1 ) THEN
344 work( k-1 ) = mid
345 ELSE
346 work( k ) = mid
347 END IF
348 i = next
349
350 100 CONTINUE
351 iter = iter + 1
352* do another loop if there are still unconverged intervals
353* However, in the last iteration, all intervals are accepted
354* since this is the best we can do.
355 IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
356*
357*
358* At this point, all the intervals have converged
359 DO 110 i = savi1, ilast
360 k = 2*i
361 ii = i - offset
362* All intervals marked by '0' have been refined.
363 IF( iwork( k-1 ).EQ.0 ) THEN
364 w( ii ) = half*( work( k-1 )+work( k ) )
365 werr( ii ) = work( k ) - w( ii )
366 END IF
367 110 CONTINUE
368*
369
370 RETURN
371*
372* End of DLARRJ
373*
374 END
dlarrj
subroutine dlarrj(n, d, e2, ifirst, ilast, rtol, offset, w, werr, work, iwork, pivmin, spdiam, info)
DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.
Definition dlarrj.f:166