1 *> \brief \b SLARRB provides limited bisection to locate eigenvalues for more accuracy.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SLARRB + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrb.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrb.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrb.f">
21 * SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
22 * RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
23 * PIVMIN, SPDIAM, TWIST, INFO )
25 * .. Scalar Arguments ..
26 * INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
27 * REAL PIVMIN, RTOL1, RTOL2, SPDIAM
29 * .. Array Arguments ..
31 * REAL D( * ), LLD( * ), W( * ),
32 * $ WERR( * ), WGAP( * ), WORK( * )
41 *> Given the relatively robust representation(RRR) L D L^T, SLARRB
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.
57 *> The order of the matrix.
62 *> D is REAL array, dimension (N)
63 *> The N diagonal elements of the diagonal matrix D.
68 *> LLD is REAL array, dimension (N-1)
69 *> The (N-1) elements L(i)*L(i)*D(i).
75 *> The index of the first eigenvalue to be computed.
81 *> The index of the last eigenvalue to be computed.
92 *> Tolerance for the convergence of the bisection intervals.
93 *> An interval [LEFT,RIGHT] has converged if
94 *> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
95 *> where GAP is the (estimated) distance to the nearest
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.
108 *> W is REAL 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 throug
112 *> On output, these estimates are refined.
115 *> \param[in,out] WGAP
117 *> WGAP is REAL 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.EQ.ILAST
121 *> then WGAP(IFIRST-OFFSET) must be set to ZERO.
122 *> On output, these gaps are refined.
125 *> \param[in,out] WERR
127 *> WERR is REAL 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.
135 *> WORK is REAL array, dimension (2*N)
141 *> IWORK is INTEGER array, dimension (2*N)
148 *> The minimum pivot in the Sturm sequence.
154 *> The spectral diameter of the matrix.
160 *> The twist index for the twisted factorization that is used
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)
176 *> \author Univ. of Tennessee
177 *> \author Univ. of California Berkeley
178 *> \author Univ. of Colorado Denver
181 *> \date September 2012
183 *> \ingroup auxOTHERauxiliary
185 *> \par Contributors:
188 *> Beresford Parlett, University of California, Berkeley, USA \n
189 *> Jim Demmel, University of California, Berkeley, USA \n
190 *> Inderjit Dhillon, University of Texas, Austin, USA \n
191 *> Osni Marques, LBNL/NERSC, USA \n
192 *> Christof Voemel, University of California, Berkeley, USA
194 * =====================================================================
195 SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
196 $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
197 $ PIVMIN, SPDIAM, TWIST, INFO )
199 * -- LAPACK auxiliary routine (version 3.4.2) --
200 * -- LAPACK is a software package provided by Univ. of Tennessee, --
201 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
204 * .. Scalar Arguments ..
205 INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
206 REAL PIVMIN, RTOL1, RTOL2, SPDIAM
208 * .. Array Arguments ..
210 REAL D( * ), LLD( * ), W( * ),
211 $ WERR( * ), WGAP( * ), WORK( * )
214 * =====================================================================
218 PARAMETER ( ZERO = 0.0E0, TWO = 2.0E0,
222 * .. Local Scalars ..
223 INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
225 REAL BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
226 $ RGAP, RIGHT, TMP, WIDTH
228 * .. External Functions ..
233 * .. Intrinsic Functions ..
234 INTRINSIC ABS, MAX, MIN
236 * .. Executable Statements ..
240 MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
242 MNWDTH = TWO * PIVMIN
245 IF((R.LT.1).OR.(R.GT.N)) R = N
247 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
248 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
249 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
250 * for an unconverged interval is set to the index of the next unconverged
251 * interval, and is -1 or 0 for a converged interval. Thus a linked
252 * list of unconverged intervals is set up.
255 * The number of unconverged intervals
257 * The last unconverged interval found
260 RGAP = WGAP( I1-OFFSET )
264 LEFT = W( II ) - WERR( II )
265 RIGHT = W( II ) + WERR( II )
268 GAP = MIN( LGAP, RGAP )
270 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
271 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
273 * Do while( NEGCNT(LEFT).GT.I-1 )
277 NEGCNT = SLANEG( N, D, LLD, LEFT, PIVMIN, R )
278 IF( NEGCNT.GT.I-1 ) THEN
284 * Do while( NEGCNT(RIGHT).LT.I )
285 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
290 NEGCNT = SLANEG( N, D, LLD, RIGHT, PIVMIN, R )
291 IF( NEGCNT.LT.I ) THEN
296 WIDTH = HALF*ABS( LEFT - RIGHT )
297 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
298 CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
299 IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
300 * This interval has already converged and does not need refinement.
301 * (Note that the gaps might change through refining the
302 * eigenvalues, however, they can only get bigger.)
303 * Remove it from the list.
305 * Make sure that I1 always points to the first unconverged interval
306 IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
307 IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
309 * unconverged interval found
320 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
321 * and while (ITER.LT.MAXITR)
329 DO 100 IP = 1, OLNINT
334 IF(II.GT.1) LGAP = WGAP( II-1 )
335 GAP = MIN( LGAP, RGAP )
339 MID = HALF*( LEFT + RIGHT )
341 * semiwidth of interval
343 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
344 CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
345 IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
346 $ ( ITER.EQ.MAXITR ) )THEN
347 * reduce number of unconverged intervals
349 * Mark interval as converged.
354 * Prev holds the last unconverged interval previously examined
355 IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
362 * Perform one bisection step
364 NEGCNT = SLANEG( N, D, LLD, MID, PIVMIN, R )
365 IF( NEGCNT.LE.I-1 ) THEN
373 * do another loop if there are still unconverged intervals
374 * However, in the last iteration, all intervals are accepted
375 * since this is the best we can do.
376 IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
379 * At this point, all the intervals have converged
380 DO 110 I = IFIRST, ILAST
383 * All intervals marked by '0' have been refined.
384 IF( IWORK( K-1 ).EQ.0 ) THEN
385 W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
386 WERR( II ) = WORK( K ) - W( II )
390 DO 111 I = IFIRST+1, ILAST
393 WGAP( II-1 ) = MAX( ZERO,
394 $ W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))