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