dlarrj()

subroutine dlarrj	(	integer	n,
		double precision, dimension( * )	d,
		double precision, dimension( * )	e2,
		integer	ifirst,
		integer	ilast,
		double precision	rtol,
		integer	offset,
		double precision, dimension( * )	w,
		double precision, dimension( * )	werr,
		double precision, dimension( * )	work,
		integer, dimension( * )	iwork,
		double precision	pivmin,
		double precision	spdiam,
		integer	info
	)

DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.

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

Purpose:

 Given the initial eigenvalue approximations of T, DLARRJ
 does  bisection to refine the eigenvalues of 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 in WERR. 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 T.
[in]	E2	E2 is DOUBLE PRECISION array, dimension (N-1) The Squares of the (N-1) subdiagonal elements of T.
[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]	RTOL	RTOL is DOUBLE PRECISION Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < RTOL*MAX(|LEFT|,|RIGHT|).
[in]	OFFSET	OFFSET is INTEGER Offset for the arrays W 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]	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 for T.
[in]	SPDIAM	SPDIAM is DOUBLE PRECISION The spectral diameter of T.
[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 165 of file dlarrj.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
      DOUBLE PRECISION   PIVMIN, RTOL, SPDIAM
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), E2( * ), W( * ),
     $                   WERR( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO, HALF
      parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
     $                   half = 0.5d0 )
      INTEGER   MAXITR
*     ..
*     .. Local Scalars ..
      INTEGER            CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
     $                   OLNINT, P, PREV, SAVI1
      DOUBLE PRECISION   DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
*
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. 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
*
*     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
      i2 = ilast
*     The number of unconverged intervals
      nint = 0
*     The last unconverged interval found
      prev = 0
      DO 75 i = i1, i2
         k = 2*i
         ii = i - offset
         left = w( ii ) - werr( ii )
         mid = w(ii)
         right = w( ii ) + werr( ii )
         width = right - mid
         tmp = max( abs( left ), abs( right ) )
 
*        The following test prevents the test of converged intervals
         IF( width.LT.rtol*tmp ) 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.i2)) i1 = i + 1
            IF((prev.GE.i1).AND.(i.LE.i2)) iwork( 2*prev-1 ) = i + 1
         ELSE
*           unconverged interval found
            prev = i
*           Make sure that [LEFT,RIGHT] contains the desired eigenvalue
*
*           Do while( CNT(LEFT).GT.I-1 )
*
            fac = one
 20         CONTINUE
            cnt = 0
            s = left
            dplus = d( 1 ) - s
            IF( dplus.LT.zero ) cnt = cnt + 1
            DO 30 j = 2, n
               dplus = d( j ) - s - e2( j-1 )/dplus
               IF( dplus.LT.zero ) cnt = cnt + 1
 30         CONTINUE
            IF( cnt.GT.i-1 ) THEN
               left = left - werr( ii )*fac
               fac = two*fac
               GO TO 20
            END IF
*
*           Do while( CNT(RIGHT).LT.I )
*
            fac = one
 50         CONTINUE
            cnt = 0
            s = right
            dplus = d( 1 ) - s
            IF( dplus.LT.zero ) cnt = cnt + 1
            DO 60 j = 2, n
               dplus = d( j ) - s - e2( j-1 )/dplus
               IF( dplus.LT.zero ) cnt = cnt + 1
 60         CONTINUE
            IF( cnt.LT.i ) THEN
               right = right + werr( ii )*fac
               fac = two*fac
               GO TO 50
            END IF
            nint = nint + 1
            iwork( k-1 ) = i + 1
            iwork( k ) = cnt
         END IF
         work( k-1 ) = left
         work( k ) = right
 75   CONTINUE
 
 
      savi1 = i1
*
*     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 p = 1, olnint
         k = 2*i
         ii = i - offset
         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 ) )
 
         IF( ( width.LT.rtol*tmp ) .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
*
         cnt = 0
         s = mid
         dplus = d( 1 ) - s
         IF( dplus.LT.zero ) cnt = cnt + 1
         DO 90 j = 2, n
            dplus = d( j ) - s - e2( j-1 )/dplus
            IF( dplus.LT.zero ) cnt = cnt + 1
 90      CONTINUE
         IF( cnt.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 = savi1, 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
*
 
      RETURN
*
*     End of DLARRJ
*

Here is the caller graph for this function: