slarrd()

subroutine slarrd	(	character	range,
		character	order,
		integer	n,
		real	vl,
		real	vu,
		integer	il,
		integer	iu,
		real, dimension( * )	gers,
		real	reltol,
		real, dimension( * )	d,
		real, dimension( * )	e,
		real, dimension( * )	e2,
		real	pivmin,
		integer	nsplit,
		integer, dimension( * )	isplit,
		integer	m,
		real, dimension( * )	w,
		real, dimension( * )	werr,
		real	wl,
		real	wu,
		integer, dimension( * )	iblock,
		integer, dimension( * )	indexw,
		real, dimension( * )	work,
		integer, dimension( * )	iwork,
		integer	info )

SLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.

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

Purpose:

!>
!> SLARRD computes the eigenvalues of a symmetric tridiagonal
!> matrix T to suitable accuracy. This is an auxiliary code to be
!> called from SSTEMR.
!> The user may ask for all eigenvalues, all eigenvalues
!> in the half-open interval (VL, VU], or the IL-th through IU-th
!> eigenvalues.
!>
!> To avoid overflow, the matrix must be scaled so that its
!> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
!> accuracy, it should not be much smaller than that.
!>
!> See W. Kahan , Report CS41, Computer Science Dept., Stanford
!> University, July 21, 1966.
!> 

Parameters

[in]	RANGE	!> RANGE is CHARACTER*1 !> = 'A': () all eigenvalues will be found. !> = 'V': () all eigenvalues in the half-open interval !> (VL, VU] will be found. !> = 'I': () the IL-th through IU-th eigenvalues (of the !> entire matrix) will be found. !>
[in]	ORDER	!> ORDER is CHARACTER*1 !> = 'B': () the eigenvalues will be grouped by !> split-off block (see IBLOCK, ISPLIT) and !> ordered from smallest to largest within !> the block. !> = 'E': () !> the eigenvalues for the entire matrix !> will be ordered from smallest to !> largest. !>
[in]	N	!> N is INTEGER !> The order of the tridiagonal matrix T. N >= 0. !>
[in]	VL	!> VL is REAL !> If RANGE='V', the lower bound of the interval to !> be searched for eigenvalues. Eigenvalues less than or equal !> to VL, or greater than VU, will not be returned. VL < VU. !> Not referenced if RANGE = 'A' or 'I'. !>
[in]	VU	!> VU is REAL !> If RANGE='V', the upper bound of the interval to !> be searched for eigenvalues. Eigenvalues less than or equal !> to VL, or greater than VU, will not be returned. VL < VU. !> Not referenced if RANGE = 'A' or 'I'. !>
[in]	IL	!> IL is INTEGER !> If RANGE='I', the index of the !> smallest eigenvalue to be returned. !> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. !> Not referenced if RANGE = 'A' or 'V'. !>
[in]	IU	!> IU is INTEGER !> If RANGE='I', the index of the !> largest eigenvalue to be returned. !> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. !> Not referenced if RANGE = 'A' or 'V'. !>
[in]	GERS	!> GERS is REAL array, dimension (2*N) !> The N Gerschgorin intervals (the i-th Gerschgorin interval !> is (GERS(2*i-1), GERS(2*i)). !>
[in]	RELTOL	!> RELTOL is REAL !> The minimum relative width of an interval. When an interval !> is narrower than RELTOL times the larger (in !> magnitude) endpoint, then it is considered to be !> sufficiently small, i.e., converged. Note: this should !> always be at least radix*machine epsilon. !>
[in]	D	!> D is REAL array, dimension (N) !> The n diagonal elements of the tridiagonal matrix T. !>
[in]	E	!> E is REAL array, dimension (N-1) !> The (n-1) off-diagonal elements of the tridiagonal matrix T. !>
[in]	E2	!> E2 is REAL array, dimension (N-1) !> The (n-1) squared off-diagonal elements of the tridiagonal matrix T. !>
[in]	PIVMIN	!> PIVMIN is REAL !> The minimum pivot allowed in the Sturm sequence for T. !>
[in]	NSPLIT	!> NSPLIT is INTEGER !> The number of diagonal blocks in the matrix T. !> 1 <= NSPLIT <= N. !>
[in]	ISPLIT	!> ISPLIT is INTEGER array, dimension (N) !> The splitting points, at which T breaks up into submatrices. !> The first submatrix consists of rows/columns 1 to ISPLIT(1), !> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), !> etc., and the NSPLIT-th consists of rows/columns !> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. !> (Only the first NSPLIT elements will actually be used, but !> since the user cannot know a priori what value NSPLIT will !> have, N words must be reserved for ISPLIT.) !>
[out]	M	!> M is INTEGER !> The actual number of eigenvalues found. 0 <= M <= N. !> (See also the description of INFO=2,3.) !>
[out]	W	!> W is REAL array, dimension (N) !> On exit, the first M elements of W will contain the !> eigenvalue approximations. SLARRD computes an interval !> I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue !> approximation is given as the interval midpoint !> W(j)= ( a_j + b_j)/2. The corresponding error is bounded by !> WERR(j) = abs( a_j - b_j)/2 !>
[out]	WERR	!> WERR is REAL array, dimension (N) !> The error bound on the corresponding eigenvalue approximation !> in W. !>
[out]	WL	!> WL is REAL !>
[out]	WU	!> WU is REAL !> The interval (WL, WU] contains all the wanted eigenvalues. !> If RANGE='V', then WL=VL and WU=VU. !> If RANGE='A', then WL and WU are the global Gerschgorin bounds !> on the spectrum. !> If RANGE='I', then WL and WU are computed by SLAEBZ from the !> index range specified. !>
[out]	IBLOCK	!> IBLOCK is INTEGER array, dimension (N) !> At each row/column j where E(j) is zero or small, the !> matrix T is considered to split into a block diagonal !> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which !> block (from 1 to the number of blocks) the eigenvalue W(i) !> belongs. (SLARRD may use the remaining N-M elements as !> workspace.) !>
[out]	INDEXW	!> INDEXW is INTEGER array, dimension (N) !> The indices of the eigenvalues within each block (submatrix); !> for example, INDEXW(i)= j and IBLOCK(i)=k imply that the !> i-th eigenvalue W(i) is the j-th eigenvalue in block k. !>
[out]	WORK	!> WORK is REAL array, dimension (4*N) !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (3*N) !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !> > 0: some or all of the eigenvalues failed to converge or !> were not computed: !> =1 or 3: Bisection failed to converge for some !> eigenvalues; these eigenvalues are flagged by a !> negative block number. The effect is that the !> eigenvalues may not be as accurate as the !> absolute and relative tolerances. This is !> generally caused by unexpectedly inaccurate !> arithmetic. !> =2 or 3: RANGE='I' only: Not all of the eigenvalues !> IL:IU were found. !> Effect: M < IU+1-IL !> Cause: non-monotonic arithmetic, causing the !> Sturm sequence to be non-monotonic. !> Cure: recalculate, using RANGE='A', and pick !> out eigenvalues IL:IU. In some cases, !> increasing the PARAMETER may !> make things work. !> = 4: RANGE='I', and the Gershgorin interval !> initially used was too small. No eigenvalues !> were computed. !> Probable cause: your machine has sloppy !> floating-point arithmetic. !> Cure: Increase the PARAMETER , !> recompile, and try again. !>

Internal Parameters:

!>  FUDGE   REAL, default = 2
!>          A  to widen the Gershgorin intervals.  Ideally,
!>          a value of 1 should work, but on machines with sloppy
!>          arithmetic, this needs to be larger.  The default for
!>          publicly released versions should be large enough to handle
!>          the worst machine around.  Note that this has no effect
!>          on accuracy of the solution.
!> 

Contributors:

W. Kahan, University of California, Berkeley, USA
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

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 323 of file slarrd.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 ..
      CHARACTER          ORDER, RANGE
      INTEGER            IL, INFO, IU, M, N, NSPLIT
      REAL                PIVMIN, RELTOL, VL, VU, WL, WU
*     ..
*     .. Array Arguments ..
      INTEGER            IBLOCK( * ), INDEXW( * ),
     $                   ISPLIT( * ), IWORK( * )
      REAL               D( * ), E( * ), E2( * ),
     $                   GERS( * ), W( * ), WERR( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO, HALF, FUDGE
      parameter( zero = 0.0e0, one = 1.0e0,
     $                     two = 2.0e0, half = one/two,
     $                     fudge = two )
      INTEGER   ALLRNG, VALRNG, INDRNG
      parameter( allrng = 1, valrng = 2, indrng = 3 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NCNVRG, TOOFEW
      INTEGER            I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
     $                   IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
     $                   ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB,
     $                   NWL, NWU
      REAL               ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2,
     $                   TNORM, UFLOW, WKILL, WLU, WUL
 
*     ..
*     .. Local Arrays ..
      INTEGER            IDUMMA( 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      REAL               SLAMCH
      EXTERNAL           lsame, ilaenv, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           slaebz
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, int, log, max, min
*     ..
*     .. Executable Statements ..
*
      info = 0
      m = 0
*
*     Quick return if possible
*
      IF( n.LE.0 ) THEN
         RETURN
      END IF
*
*     Decode RANGE
*
      IF( lsame( range, 'A' ) ) THEN
         irange = allrng
      ELSE IF( lsame( range, 'V' ) ) THEN
         irange = valrng
      ELSE IF( lsame( range, 'I' ) ) THEN
         irange = indrng
      ELSE
         irange = 0
      END IF
*
*     Check for Errors
*
      IF( irange.LE.0 ) THEN
         info = -1
      ELSE IF( .NOT.(lsame(order,'B').OR.lsame(order,'E')) ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( irange.EQ.valrng ) THEN
         IF( vl.GE.vu )
     $      info = -5
      ELSE IF( irange.EQ.indrng .AND.
     $        ( il.LT.1 .OR. il.GT.max( 1, n ) ) ) THEN
         info = -6
      ELSE IF( irange.EQ.indrng .AND.
     $        ( iu.LT.min( n, il ) .OR. iu.GT.n ) ) THEN
         info = -7
      END IF
*
      IF( info.NE.0 ) THEN
         RETURN
      END IF
 
*     Initialize error flags
      ncnvrg = .false.
      toofew = .false.
 
*     Simplification:
      IF( irange.EQ.indrng .AND. il.EQ.1 .AND. iu.EQ.n ) irange = 1
 
*     Get machine constants
      eps = slamch( 'P' )
      uflow = slamch( 'U' )
 
 
*     Special Case when N=1
*     Treat case of 1x1 matrix for quick return
      IF( n.EQ.1 ) THEN
         IF( (irange.EQ.allrng).OR.
     $       ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
     $       ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
            m = 1
            w(1) = d(1)
*           The computation error of the eigenvalue is zero
            werr(1) = zero
            iblock( 1 ) = 1
            indexw( 1 ) = 1
         ENDIF
         RETURN
      END IF
 
*     NB is the minimum vector length for vector bisection, or 0
*     if only scalar is to be done.
      nb = ilaenv( 1, 'SSTEBZ', ' ', n, -1, -1, -1 )
      IF( nb.LE.1 ) nb = 0
 
*     Find global spectral radius
      gl = d(1)
      gu = d(1)
      DO 5 i = 1,n
         gl =  min( gl, gers( 2*i - 1))
         gu = max( gu, gers(2*i) )
 5    CONTINUE
*     Compute global Gerschgorin bounds and spectral diameter
      tnorm = max( abs( gl ), abs( gu ) )
      gl = gl - fudge*tnorm*eps*real( n ) - fudge*two*pivmin
      gu = gu + fudge*tnorm*eps*real( n ) + fudge*two*pivmin
*     [JAN/28/2009] remove the line below since SPDIAM variable not use
*     SPDIAM = GU - GL
*     Input arguments for SLAEBZ:
*     The relative tolerance.  An interval (a,b] lies within
*     "relative tolerance" if  b-a < RELTOL*max(|a|,|b|),
      rtoli = reltol
*     Set the absolute tolerance for interval convergence to zero to force
*     interval convergence based on relative size of the interval.
*     This is dangerous because intervals might not converge when RELTOL is
*     small. But at least a very small number should be selected so that for
*     strongly graded matrices, the code can get relatively accurate
*     eigenvalues.
      atoli = fudge*two*uflow + fudge*two*pivmin
 
      IF( irange.EQ.indrng ) THEN
 
*        RANGE='I': Compute an interval containing eigenvalues
*        IL through IU. The initial interval [GL,GU] from the global
*        Gerschgorin bounds GL and GU is refined by SLAEBZ.
         itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
     $           log( two ) ) + 2
         work( n+1 ) = gl
         work( n+2 ) = gl
         work( n+3 ) = gu
         work( n+4 ) = gu
         work( n+5 ) = gl
         work( n+6 ) = gu
         iwork( 1 ) = -1
         iwork( 2 ) = -1
         iwork( 3 ) = n + 1
         iwork( 4 ) = n + 1
         iwork( 5 ) = il - 1
         iwork( 6 ) = iu
*
         CALL slaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin,
     $         d, e, e2, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
     $                iwork, w, iblock, iinfo )
         IF( iinfo .NE. 0 ) THEN
            info = iinfo
            RETURN
         END IF
*        On exit, output intervals may not be ordered by ascending negcount
         IF( iwork( 6 ).EQ.iu ) THEN
            wl = work( n+1 )
            wlu = work( n+3 )
            nwl = iwork( 1 )
            wu = work( n+4 )
            wul = work( n+2 )
            nwu = iwork( 4 )
         ELSE
            wl = work( n+2 )
            wlu = work( n+4 )
            nwl = iwork( 2 )
            wu = work( n+3 )
            wul = work( n+1 )
            nwu = iwork( 3 )
         END IF
*        On exit, the interval [WL, WLU] contains a value with negcount NWL,
*        and [WUL, WU] contains a value with negcount NWU.
         IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
            info = 4
            RETURN
         END IF
 
      ELSEIF( irange.EQ.valrng ) THEN
         wl = vl
         wu = vu
 
      ELSEIF( irange.EQ.allrng ) THEN
         wl = gl
         wu = gu
      ENDIF
 
 
 
*     Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
*     NWL accumulates the number of eigenvalues .le. WL,
*     NWU accumulates the number of eigenvalues .le. WU
      m = 0
      iend = 0
      info = 0
      nwl = 0
      nwu = 0
*
      DO 70 jblk = 1, nsplit
         ioff = iend
         ibegin = ioff + 1
         iend = isplit( jblk )
         in = iend - ioff
*
         IF( in.EQ.1 ) THEN
*           1x1 block
            IF( wl.GE.d( ibegin )-pivmin )
     $         nwl = nwl + 1
            IF( wu.GE.d( ibegin )-pivmin )
     $         nwu = nwu + 1
            IF( irange.EQ.allrng .OR.
     $           ( wl.LT.d( ibegin )-pivmin
     $             .AND. wu.GE. d( ibegin )-pivmin ) ) THEN
               m = m + 1
               w( m ) = d( ibegin )
               werr(m) = zero
*              The gap for a single block doesn't matter for the later
*              algorithm and is assigned an arbitrary large value
               iblock( m ) = jblk
               indexw( m ) = 1
            END IF
 
*        Disabled 2x2 case because of a failure on the following matrix
*        RANGE = 'I', IL = IU = 4
*          Original Tridiagonal, d = [
*           -0.150102010615740E+00
*           -0.849897989384260E+00
*           -0.128208148052635E-15
*            0.128257718286320E-15
*          ];
*          e = [
*           -0.357171383266986E+00
*           -0.180411241501588E-15
*           -0.175152352710251E-15
*          ];
*
*         ELSE IF( IN.EQ.2 ) THEN
**           2x2 block
*            DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
*            TMP1 = HALF*(D(IBEGIN)+D(IEND))
*            L1 = TMP1 - DISC
*            IF( WL.GE. L1-PIVMIN )
*     $         NWL = NWL + 1
*            IF( WU.GE. L1-PIVMIN )
*     $         NWU = NWU + 1
*            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
*     $          L1-PIVMIN ) ) THEN
*               M = M + 1
*               W( M ) = L1
**              The uncertainty of eigenvalues of a 2x2 matrix is very small
*               WERR( M ) = EPS * ABS( W( M ) ) * TWO
*               IBLOCK( M ) = JBLK
*               INDEXW( M ) = 1
*            ENDIF
*            L2 = TMP1 + DISC
*            IF( WL.GE. L2-PIVMIN )
*     $         NWL = NWL + 1
*            IF( WU.GE. L2-PIVMIN )
*     $         NWU = NWU + 1
*            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
*     $          L2-PIVMIN ) ) THEN
*               M = M + 1
*               W( M ) = L2
**              The uncertainty of eigenvalues of a 2x2 matrix is very small
*               WERR( M ) = EPS * ABS( W( M ) ) * TWO
*               IBLOCK( M ) = JBLK
*               INDEXW( M ) = 2
*            ENDIF
         ELSE
*           General Case - block of size IN >= 2
*           Compute local Gerschgorin interval and use it as the initial
*           interval for SLAEBZ
            gu = d( ibegin )
            gl = d( ibegin )
            tmp1 = zero
 
            DO 40 j = ibegin, iend
               gl =  min( gl, gers( 2*j - 1))
               gu = max( gu, gers(2*j) )
   40       CONTINUE
*           [JAN/28/2009]
*           change SPDIAM by TNORM in lines 2 and 3 thereafter
*           line 1: remove computation of SPDIAM (not useful anymore)
*           SPDIAM = GU - GL
*           GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
*           GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
            gl = gl - fudge*tnorm*eps*real( in ) - fudge*pivmin
            gu = gu + fudge*tnorm*eps*real( in ) + fudge*pivmin
*
            IF( irange.GT.1 ) THEN
               IF( gu.LT.wl ) THEN
*                 the local block contains none of the wanted eigenvalues
                  nwl = nwl + in
                  nwu = nwu + in
                  GO TO 70
               END IF
*              refine search interval if possible, only range (WL,WU] matters
               gl = max( gl, wl )
               gu = min( gu, wu )
               IF( gl.GE.gu )
     $            GO TO 70
            END IF
 
*           Find negcount of initial interval boundaries GL and GU
            work( n+1 ) = gl
            work( n+in+1 ) = gu
            CALL slaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
     $                   d( ibegin ), e( ibegin ), e2( ibegin ),
     $                   idumma, work( n+1 ), work( n+2*in+1 ), im,
     $                   iwork, w( m+1 ), iblock( m+1 ), iinfo )
            IF( iinfo .NE. 0 ) THEN
               info = iinfo
               RETURN
            END IF
*
            nwl = nwl + iwork( 1 )
            nwu = nwu + iwork( in+1 )
            iwoff = m - iwork( 1 )
 
*           Compute Eigenvalues
            itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
     $              log( two ) ) + 2
            CALL slaebz( 2, itmax, in, in, 1, nb, atoli, rtoli,
     $                   pivmin,
     $                   d( ibegin ), e( ibegin ), e2( ibegin ),
     $                   idumma, work( n+1 ), work( n+2*in+1 ), iout,
     $                   iwork, w( m+1 ), iblock( m+1 ), iinfo )
            IF( iinfo .NE. 0 ) THEN
               info = iinfo
               RETURN
            END IF
*
*           Copy eigenvalues into W and IBLOCK
*           Use -JBLK for block number for unconverged eigenvalues.
*           Loop over the number of output intervals from SLAEBZ
            DO 60 j = 1, iout
*              eigenvalue approximation is middle point of interval
               tmp1 = half*( work( j+n )+work( j+in+n ) )
*              semi length of error interval
               tmp2 = half*abs( work( j+n )-work( j+in+n ) )
               IF( j.GT.iout-iinfo ) THEN
*                 Flag non-convergence.
                  ncnvrg = .true.
                  ib = -jblk
               ELSE
                  ib = jblk
               END IF
               DO 50 je = iwork( j ) + 1 + iwoff,
     $                 iwork( j+in ) + iwoff
                  w( je ) = tmp1
                  werr( je ) = tmp2
                  indexw( je ) = je - iwoff
                  iblock( je ) = ib
   50          CONTINUE
   60       CONTINUE
*
            m = m + im
         END IF
   70 CONTINUE
 
*     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
*     If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
      IF( irange.EQ.indrng ) THEN
         idiscl = il - 1 - nwl
         idiscu = nwu - iu
*
         IF( idiscl.GT.0 ) THEN
            im = 0
            DO 80 je = 1, m
*              Remove some of the smallest eigenvalues from the left so that
*              at the end IDISCL =0. Move all eigenvalues up to the left.
               IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
                  idiscl = idiscl - 1
               ELSE
                  im = im + 1
                  w( im ) = w( je )
                  werr( im ) = werr( je )
                  indexw( im ) = indexw( je )
                  iblock( im ) = iblock( je )
               END IF
 80         CONTINUE
            m = im
         END IF
         IF( idiscu.GT.0 ) THEN
*           Remove some of the largest eigenvalues from the right so that
*           at the end IDISCU =0. Move all eigenvalues up to the left.
            im=m+1
            DO 81 je = m, 1, -1
               IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
                  idiscu = idiscu - 1
               ELSE
                  im = im - 1
                  w( im ) = w( je )
                  werr( im ) = werr( je )
                  indexw( im ) = indexw( je )
                  iblock( im ) = iblock( je )
               END IF
 81         CONTINUE
            jee = 0
            DO 82 je = im, m
               jee = jee + 1
               w( jee ) = w( je )
               werr( jee ) = werr( je )
               indexw( jee ) = indexw( je )
               iblock( jee ) = iblock( je )
 82         CONTINUE
            m = m-im+1
         END IF
 
         IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
*           Code to deal with effects of bad arithmetic. (If N(w) is
*           monotone non-decreasing, this should never happen.)
*           Some low eigenvalues to be discarded are not in (WL,WLU],
*           or high eigenvalues to be discarded are not in (WUL,WU]
*           so just kill off the smallest IDISCL/largest IDISCU
*           eigenvalues, by marking the corresponding IBLOCK = 0
            IF( idiscl.GT.0 ) THEN
               wkill = wu
               DO 100 jdisc = 1, idiscl
                  iw = 0
                  DO 90 je = 1, m
                     IF( iblock( je ).NE.0 .AND.
     $                    ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
                        iw = je
                        wkill = w( je )
                     END IF
 90               CONTINUE
                  iblock( iw ) = 0
 100           CONTINUE
            END IF
            IF( idiscu.GT.0 ) THEN
               wkill = wl
               DO 120 jdisc = 1, idiscu
                  iw = 0
                  DO 110 je = 1, m
                     IF( iblock( je ).NE.0 .AND.
     $                    ( w( je ).GE.wkill .OR. iw.EQ.0 ) ) THEN
                        iw = je
                        wkill = w( je )
                     END IF
 110              CONTINUE
                  iblock( iw ) = 0
 120           CONTINUE
            END IF
*           Now erase all eigenvalues with IBLOCK set to zero
            im = 0
            DO 130 je = 1, m
               IF( iblock( je ).NE.0 ) THEN
                  im = im + 1
                  w( im ) = w( je )
                  werr( im ) = werr( je )
                  indexw( im ) = indexw( je )
                  iblock( im ) = iblock( je )
               END IF
 130        CONTINUE
            m = im
         END IF
         IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
            toofew = .true.
         END IF
      END IF
*
      IF(( irange.EQ.allrng .AND. m.NE.n ).OR.
     $   ( irange.EQ.indrng .AND. m.NE.iu-il+1 ) ) THEN
         toofew = .true.
      END IF
 
*     If ORDER='B', do nothing the eigenvalues are already sorted by
*        block.
*     If ORDER='E', sort the eigenvalues from smallest to largest
 
      IF( lsame(order,'E') .AND. nsplit.GT.1 ) THEN
         DO 150 je = 1, m - 1
            ie = 0
            tmp1 = w( je )
            DO 140 j = je + 1, m
               IF( w( j ).LT.tmp1 ) THEN
                  ie = j
                  tmp1 = w( j )
               END IF
  140       CONTINUE
            IF( ie.NE.0 ) THEN
               tmp2 = werr( ie )
               itmp1 = iblock( ie )
               itmp2 = indexw( ie )
               w( ie ) = w( je )
               werr( ie ) = werr( je )
               iblock( ie ) = iblock( je )
               indexw( ie ) = indexw( je )
               w( je ) = tmp1
               werr( je ) = tmp2
               iblock( je ) = itmp1
               indexw( je ) = itmp2
            END IF
  150    CONTINUE
      END IF
*
      info = 0
      IF( ncnvrg )
     $   info = info + 1
      IF( toofew )
     $   info = info + 2
      RETURN
*
*     End of SLARRD
*

Here is the call graph for this function:

Here is the caller graph for this function: