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