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.

Download DLARRB + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> Given the relatively robust representation(RRR) L D L^T, DLARRB
!> does  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. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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 191 of file dlarrb.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            IFIRST, ILAST, INFO, N, OFFSET, TWIST
      DOUBLE PRECISION   PIVMIN, RTOL1, RTOL2, SPDIAM
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), LLD( * ), W( * ),
     $                   WERR( * ), WGAP( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, TWO, HALF
      parameter( zero = 0.0d0, two = 2.0d0,
     $                   half = 0.5d0 )
      INTEGER   MAXITR
*     ..
*     .. Local Scalars ..
      INTEGER            I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
     $                   OLNINT, PREV, R
      DOUBLE PRECISION   BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
     $                   RGAP, RIGHT, TMP, WIDTH
*     ..
*     .. External Functions ..
      INTEGER            DLANEG
      EXTERNAL           dlaneg
*
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Quick return if possible
*
      IF( n.LE.0 ) THEN
         RETURN
      END IF
*
      maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
     $           log( two ) ) + 2
      mnwdth = two * pivmin
*
      r = twist
      IF((r.LT.1).OR.(r.GT.n)) r = n
*
*     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
*     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
*     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
*     for an unconverged interval is set to the index of the next unconverged
*     interval, and is -1 or 0 for a converged interval. Thus a linked
*     list of unconverged intervals is set up.
*
      i1 = ifirst
*     The number of unconverged intervals
      nint = 0
*     The last unconverged interval found
      prev = 0
 
      rgap = wgap( i1-offset )
      DO 75 i = i1, ilast
         k = 2*i
         ii = i - offset
         left = w( ii ) - werr( ii )
         right = w( ii ) + werr( ii )
         lgap = rgap
         rgap = wgap( ii )
         gap = min( lgap, rgap )
 
*        Make sure that [LEFT,RIGHT] contains the desired eigenvalue
*        Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
*
*        Do while( NEGCNT(LEFT).GT.I-1 )
*
         back = werr( ii )
 20      CONTINUE
         negcnt = dlaneg( n, d, lld, left, pivmin, r )
         IF( negcnt.GT.i-1 ) THEN
            left = left - back
            back = two*back
            GO TO 20
         END IF
*
*        Do while( NEGCNT(RIGHT).LT.I )
*        Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
*
         back = werr( ii )
 50      CONTINUE
 
         negcnt = dlaneg( n, d, lld, right, pivmin, r )
          IF( negcnt.LT.i ) THEN
             right = right + back
             back = two*back
             GO TO 50
          END IF
         width = half*abs( left - right )
         tmp = max( abs( left ), abs( right ) )
         cvrgd = max(rtol1*gap,rtol2*tmp)
         IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
*           This interval has already converged and does not need refinement.
*           (Note that the gaps might change through refining the
*            eigenvalues, however, they can only get bigger.)
*           Remove it from the list.
            iwork( k-1 ) = -1
*           Make sure that I1 always points to the first unconverged interval
            IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
            IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
         ELSE
*           unconverged interval found
            prev = i
            nint = nint + 1
            iwork( k-1 ) = i + 1
            iwork( k ) = negcnt
         END IF
         work( k-1 ) = left
         work( k ) = right
 75   CONTINUE
 
*
*     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
*     and while (ITER.LT.MAXITR)
*
      iter = 0
 80   CONTINUE
      prev = i1 - 1
      i = i1
      olnint = nint
 
      DO 100 ip = 1, olnint
         k = 2*i
         ii = i - offset
         rgap = wgap( ii )
         lgap = rgap
         IF(ii.GT.1) lgap = wgap( ii-1 )
         gap = min( lgap, rgap )
         next = iwork( k-1 )
         left = work( k-1 )
         right = work( k )
         mid = half*( left + right )
 
*        semiwidth of interval
         width = right - mid
         tmp = max( abs( left ), abs( right ) )
         cvrgd = max(rtol1*gap,rtol2*tmp)
         IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
     $       ( iter.EQ.maxitr ) )THEN
*           reduce number of unconverged intervals
            nint = nint - 1
*           Mark interval as converged.
            iwork( k-1 ) = 0
            IF( i1.EQ.i ) THEN
               i1 = next
            ELSE
*              Prev holds the last unconverged interval previously examined
               IF(prev.GE.i1) iwork( 2*prev-1 ) = next
            END IF
            i = next
            GO TO 100
         END IF
         prev = i
*
*        Perform one bisection step
*
         negcnt = dlaneg( n, d, lld, mid, pivmin, r )
         IF( negcnt.LE.i-1 ) THEN
            work( k-1 ) = mid
         ELSE
            work( k ) = mid
         END IF
         i = next
 100  CONTINUE
      iter = iter + 1
*     do another loop if there are still unconverged intervals
*     However, in the last iteration, all intervals are accepted
*     since this is the best we can do.
      IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
*
*
*     At this point, all the intervals have converged
      DO 110 i = ifirst, ilast
         k = 2*i
         ii = i - offset
*        All intervals marked by '0' have been refined.
         IF( iwork( k-1 ).EQ.0 ) THEN
            w( ii ) = half*( work( k-1 )+work( k ) )
            werr( ii ) = work( k ) - w( ii )
         END IF
 110  CONTINUE
*
      DO 111 i = ifirst+1, ilast
         k = 2*i
         ii = i - offset
         wgap( ii-1 ) = max( zero,
     $                     w(ii) - werr(ii) - w( ii-1 ) - werr( ii-1 ))
 111  CONTINUE
 
      RETURN
*
*     End of DLARRB
*

Here is the caller graph for this function: