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