slarrj()

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

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

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

Purpose:

!>
!> Given the initial eigenvalue approximations of T, SLARRJ
!> 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 REAL array, dimension (N) !> The N diagonal elements of T. !>
[in]	E2	!> E2 is REAL 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 REAL !> 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 REAL 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 REAL 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 REAL array, dimension (2*N) !> Workspace. !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (2*N) !> Workspace. !>
[in]	PIVMIN	!> PIVMIN is REAL !> The minimum pivot in the Sturm sequence for T. !>
[in]	SPDIAM	!> SPDIAM is REAL !> 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 163 of file slarrj.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
      REAL               PIVMIN, RTOL, SPDIAM
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      REAL               D( * ), E2( * ), W( * ),
     $                   WERR( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO, HALF
      parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
     $                   half = 0.5e0 )
      INTEGER   MAXITR
*     ..
*     .. Local Scalars ..
      INTEGER            CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
     $                   OLNINT, P, PREV, SAVI1
      REAL               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 SLARRJ
*

Here is the caller graph for this function: