dlarrb2()

subroutine dlarrb2	(	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	lgpvmn,
		double precision	lgspdm,
		integer	twist,
		integer	info
	)

Definition at line 1 of file dlarrb2.f.

*
*  -- ScaLAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, Univ of Colorado Denver
*     July 4, 2010
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      INTEGER            IFIRST, ILAST, INFO, N, OFFSET, TWIST
      DOUBLE PRECISION   LGPVMN, LGSPDM, PIVMIN, 
     $                   RTOL1, RTOL2
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), LLD( * ), W( * ),
     $                   WERR( * ), WGAP( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  Given the relatively robust representation(RRR) L D L^T, DLARRB2
*  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.
*
*  NOTE: 
*  There are very few minor differences between DLARRB from LAPACK
*  and this current subroutine DLARRB2.
*  The most important reason for creating this nearly identical copy
*  is profiling: in the ScaLAPACK MRRR algorithm, eigenvalue computation 
*  using DLARRB2 is used for refinement in the construction of 
*  the representation tree, as opposed to the initial computation of the
*  eigenvalues for the root RRR which uses DLARRB. When profiling,
*  this allows an easy quantification of refinement work vs. computing
*  eigenvalues of the root.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix.
*
*  D       (input) DOUBLE PRECISION array, dimension (N)
*          The N diagonal elements of the diagonal matrix D.
*
*  LLD     (input) DOUBLE PRECISION array, dimension (N-1)
*          The (N-1) elements L(i)*L(i)*D(i).
*
*  IFIRST  (input) INTEGER
*          The index of the first eigenvalue to be computed.
*
*  ILAST   (input) INTEGER
*          The index of the last eigenvalue to be computed.
*
*  RTOL1   (input) DOUBLE PRECISION
*  RTOL2   (input) DOUBLE PRECISION
*          Tolerance for the convergence of the bisection intervals.
*          An interval [LEFT,RIGHT] has converged if
*          RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*          where GAP is the (estimated) distance to the nearest
*          eigenvalue.
*
*  OFFSET  (input) 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.
*
*  W       (input/output) 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.
*
*  WGAP    (input/output) 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.EQ.ILAST
*          then WGAP(IFIRST-OFFSET) must be set to ZERO.
*          On output, these gaps are refined.
*
*  WERR    (input/output) 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.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
*          Workspace.
*
*  IWORK   (workspace) INTEGER array, dimension (2*N)
*          Workspace.
*
*  PIVMIN  (input) DOUBLE PRECISION 
*          The minimum pivot in the sturm sequence.
*
*  LGPVMN  (input) DOUBLE PRECISION
*          Logarithm of PIVMIN, precomputed.
*
*  LGSPDM  (input) DOUBLE PRECISION 
*          Logarithm of the spectral diameter, precomputed.
*
*  TWIST   (input) 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)
*
*  INFO    (output) INTEGER
*          Error flag.
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, TWO, HALF
      parameter( zero = 0.0d0, two = 2.0d0,
     $                   half = 0.5d0 )
      INTEGER   MAXITR
*     ..
*     .. Local Scalars ..
      INTEGER            I, I1, II, INDLLD, IP, ITER, J, K, NEGCNT,
     $                   NEXT, NINT, OLNINT, PREV, R
      DOUBLE PRECISION   BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
     $                   RGAP, RIGHT, SAVGAP, TMP, WIDTH
      LOGICAL   PARANOID
*     ..
*     .. External Functions ..
      LOGICAL            DISNAN
      DOUBLE PRECISION   DLAMCH
      INTEGER            DLANEG2A
      EXTERNAL           disnan, dlamch, 
     $                   dlaneg2a
*
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*     ..
*     .. Executable Statements ..
*
      info = 0
*     
*     Turn on paranoid check for rounding errors 
*     invalidating uncertainty intervals of eigenvalues
*
      paranoid = .true.
*
      maxitr = int( ( lgspdm - lgpvmn ) / log( two ) ) + 2
      mnwdth = two * pivmin
*
      r = twist
*
      indlld = 2*n     
      DO 5 j = 1, n-1 
         i=2*j
         work(indlld+i-1) = d(j)
         work(indlld+i) = lld(j)
  5   CONTINUE
      work(indlld+2*n-1) = d(n)
*
      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 )
 
         IF((abs(left).LE.16*pivmin).OR.(abs(right).LE.16*pivmin))
     $      THEN
            info = -1
            RETURN
         ENDIF
 
         IF( paranoid ) THEN
*        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 = dlaneg2a( n, work(indlld+1), 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 = dlaneg2a( n, work(indlld+1),right, pivmin, r )
 
         IF( negcnt.LT.i ) THEN
             right = right + back
             back = two*back
             GO TO 50
         END IF
         ENDIF
 
         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 = dlaneg2a( n, work(indlld+1), 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
*
*     save this gap to restore it after the loop
      savgap = wgap( ilast-offset )
*
      left = work( 2*ifirst-1 )
      DO 110 i = ifirst, ilast
         k = 2*i
         ii = i - offset
*        RIGHT is the right boundary of this current interval
         right = work( k ) 
*        All intervals marked by '0' have been refined.
         IF( iwork( k-1 ).EQ.0 ) THEN
            w( ii ) = half*( left+right )
            werr( ii ) = right - w( ii )
         END IF
*        Left is the boundary of the next interval
         left = work( k +1 ) 
         wgap( ii ) = max( zero, left - right )
 110  CONTINUE
*     restore the last gap which was overwritten by garbage
      wgap( ilast-offset ) = savgap
 
      RETURN
*
*     End of DLARRB2
*

Here is the caller graph for this function: