 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dlarrb()

 subroutine dlarrb ( integer N, double precision, dimension( * ) D, double precision, dimension( * ) LLD, integer IFIRST, integer ILAST, double precision RTOL1, double precision RTOL2, integer OFFSET, double precision, dimension( * ) W, double precision, dimension( * ) WGAP, double precision, dimension( * ) WERR, double precision, dimension( * ) WORK, integer, dimension( * ) IWORK, double precision PIVMIN, double precision SPDIAM, integer TWIST, integer INFO )

DLARRB provides limited bisection to locate eigenvalues for more accuracy.

Purpose:
``` Given the relatively robust representation(RRR) L D L^T, DLARRB
does "limited" bisection to refine the eigenvalues of L D L^T,
W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
guesses for these eigenvalues are input in W, the corresponding estimate
of the error in these guesses and their gaps are input in WERR
and WGAP, respectively. During bisection, intervals
[left, right] are maintained by storing their mid-points and
semi-widths in the arrays W and WERR respectively.```
Parameters
 [in] N ``` N is INTEGER The order of the matrix.``` [in] D ``` D is DOUBLE PRECISION array, dimension (N) The N diagonal elements of the diagonal matrix D.``` [in] LLD ``` LLD is DOUBLE PRECISION array, dimension (N-1) The (N-1) elements L(i)*L(i)*D(i).``` [in] IFIRST ``` IFIRST is INTEGER The index of the first eigenvalue to be computed.``` [in] ILAST ``` ILAST is INTEGER The index of the last eigenvalue to be computed.``` [in] RTOL1 ` RTOL1 is DOUBLE PRECISION` [in] RTOL2 ``` RTOL2 is DOUBLE PRECISION Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) where GAP is the (estimated) distance to the nearest eigenvalue.``` [in] OFFSET ``` OFFSET is INTEGER Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET through ILAST-OFFSET elements of these arrays are to be used.``` [in,out] W ``` W is DOUBLE PRECISION array, dimension (N) On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST. On output, these estimates are refined.``` [in,out] WGAP ``` WGAP is DOUBLE PRECISION array, dimension (N-1) On input, the (estimated) gaps between consecutive eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between eigenvalues I and I+1. Note that if IFIRST = ILAST then WGAP(IFIRST-OFFSET) must be set to ZERO. On output, these gaps are refined.``` [in,out] WERR ``` WERR is DOUBLE PRECISION array, dimension (N) On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are the errors in the estimates of the corresponding elements in W. On output, these errors are refined.``` [out] WORK ``` WORK is DOUBLE PRECISION array, dimension (2*N) Workspace.``` [out] IWORK ``` IWORK is INTEGER array, dimension (2*N) Workspace.``` [in] PIVMIN ``` PIVMIN is DOUBLE PRECISION The minimum pivot in the Sturm sequence.``` [in] SPDIAM ``` SPDIAM is DOUBLE PRECISION The spectral diameter of the matrix.``` [in] TWIST ``` TWIST is INTEGER The twist index for the twisted factorization that is used for the negcount. TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)``` [out] INFO ``` INFO is INTEGER Error flag.```
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 193 of file dlarrb.f.

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 *
integer function dlaneg(N, D, LLD, SIGMA, PIVMIN, R)
DLANEG computes the Sturm count.
Definition: dlaneg.f:118
Here is the caller graph for this function: