slarrd2.f

Go to the documentation of this file.

      SUBROUTINE slarrd2( RANGE, ORDER, N, VL, VU, IL, IU, GERS, 
     $                    RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, 
     $                    M, W, WERR, WL, WU, IBLOCK, INDEXW, 
     $                    WORK, IWORK, DOL, DOU, INFO )
*
*  -- ScaLAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
*     July 4, 2010
*
*     .. Scalar Arguments ..
      CHARACTER          ORDER, RANGE
      INTEGER            DOL, DOU, 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( * )
*     ..
*
*  Purpose
*  =======
*
*  SLARRD2 computes the eigenvalues of a symmetric tridiagonal
*  matrix T to limited initial accuracy. This is an auxiliary code to be 
*  called from SLARRE2A.
* 
*  SLARRD2 has been created using the LAPACK code SLARRD
*  which itself stems from SSTEBZ. The motivation for creating
*  SLARRD2 is efficiency: When computing eigenvalues in parallel
*  and the input tridiagonal matrix splits into blocks, SLARRD2 
*  can skip over blocks which contain none of the eigenvalues from
*  DOL to DOU for which the processor responsible. In extreme cases (such
*  as large matrices consisting of many blocks of small size, e.g. 2x2,
*  the gain can be substantial.
*
*  Arguments
*  =========
*
*  RANGE   (input) CHARACTER
*          = 'A': ("All")   all eigenvalues will be found.
*          = 'V': ("Value") all eigenvalues in the half-open interval
*                           (VL, VU] will be found.
*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
*                           entire matrix) will be found.
*
*  ORDER   (input) CHARACTER
*          = 'B': ("By Block") the eigenvalues will be grouped by
*                              split-off block (see IBLOCK, ISPLIT) and
*                              ordered from smallest to largest within
*                              the block.
*          = 'E': ("Entire matrix")
*                              the eigenvalues for the entire matrix
*                              will be ordered from smallest to
*                              largest.
*
*  N       (input) INTEGER
*          The order of the tridiagonal matrix T.  N >= 0.
*
*  VL      (input) REAL            
*  VU      (input) REAL            
*          If RANGE='V', the lower and upper bounds 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'.
*
*  IL      (input) INTEGER
*  IU      (input) INTEGER
*          If RANGE='I', the indices (in ascending order) of the
*          smallest and largest eigenvalues 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'.
*
*  GERS    (input) REAL             array, dimension (2*N)
*          The N Gerschgorin intervals (the i-th Gerschgorin interval
*          is (GERS(2*i-1), GERS(2*i)).
*
*  RELTOL  (input) 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.
*
*  D       (input) REAL             array, dimension (N)
*          The n diagonal elements of the tridiagonal matrix T.
*
*  E       (input) REAL             array, dimension (N-1)
*          The (n-1) off-diagonal elements of the tridiagonal matrix T.
*
*  E2      (input) REAL             array, dimension (N-1)
*          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
*
*  PIVMIN  (input) REAL            
*          The minimum pivot allowed in the sturm sequence for T.
*
*  NSPLIT  (input) INTEGER
*          The number of diagonal blocks in the matrix T.
*          1 <= NSPLIT <= N.
*
*  ISPLIT  (input) 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.)
*
*  M       (output) INTEGER
*          The actual number of eigenvalues found. 0 <= M <= N.
*          (See also the description of INFO=2,3.)
*
*  W       (output) REAL             array, dimension (N)
*          On exit, the first M elements of W will contain the
*          eigenvalue approximations. SLARRD2 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      
*
*  WERR    (output) REAL             array, dimension (N)
*          The error bound on the corresponding eigenvalue approximation 
*          in W.
*
*  WL      (output) REAL            
*  WU      (output) 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.
*                        
*  IBLOCK  (output) 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.  (SLARRD2 may use the remaining N-M elements as
*          workspace.)
*
*  INDEXW  (output) 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.
*
*  WORK    (workspace) REAL             array, dimension (4*N)
*
*  IWORK   (workspace) INTEGER array, dimension (3*N)
*
*  DOL     (input) INTEGER
*  DOU     (input) INTEGER
*          If the user wants to work on only a selected part of the 
*          representation tree, he can specify an index range DOL:DOU.
*          Otherwise, the setting DOL=1, DOU=N should be applied. 
*          Note that DOL and DOU refer to the order in which the eigenvalues 
*          are stored in W. 
*
*  INFO    (output) 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 "FUDGE" 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 "FUDGE",
*                              recompile, and try again.
*
*  Internal Parameters
*  ===================
*
*  FUDGE   REAL            , default = 2 originally, increased to 10.
*          A "fudge factor" 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.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO, HALF, FUDGE
      PARAMETER          ( ZERO = 0.0e0, one = 1.0e0, 
     $                     two = 2.0e0, half = one/two,
     $                     fudge = 10.0e0 )
 
*     ..
*     .. Local Scalars ..
      LOGICAL            NCNVRG, TOOFEW
      INTEGER            I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
     $                   IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
     $                   itmp1, itmp2, iw, iwoff, j, jblk, jdisc, je,
     $                   jee, nb, nwl, nwu
      REAL               ATOLI, EPS, GL, GU, RTOLI, SPDIAM, 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, sqrt
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Decode RANGE
*
      IF( lsame( range, 'A' ) ) THEN
         irange = 1
      ELSE IF( lsame( range, 'V' ) ) THEN
         irange = 2
      ELSE IF( lsame( range, 'i' ) ) THEN
         IRANGE = 3
      ELSE
         IRANGE = 0
      END IF
*
*     Decode ORDER
*
      IF( LSAME( ORDER, 'b' ) ) THEN
         IORDER = 2
      ELSE IF( LSAME( ORDER, 'e' ) ) THEN
         IORDER = 1
      ELSE
         IORDER = 0
      END IF
*
*     Check for Errors
*
.LE.      IF( IRANGE0 ) THEN
         INFO = -1
.LE.      ELSE IF( IORDER0 ) THEN
         INFO = -2
.LT.      ELSE IF( N0 ) THEN
         INFO = -3
.EQ.      ELSE IF( IRANGE2 ) THEN
.GE.         IF( VLVU )
     $      INFO = -5
.EQ..AND..LT..OR..GT.      ELSE IF( IRANGE3  ( IL1  ILMAX( 1, N ) ) )
     $          THEN
         INFO = -6
.EQ..AND..LT..OR..GT.      ELSE IF( IRANGE3  ( IUMIN( N, IL )  IUN ) )
     $          THEN
         INFO = -7
      END IF
*
.NE.      IF( INFO0 ) THEN
         RETURN
      END IF
 
*     Initialize error flags
      INFO = 0
      NCNVRG = .FALSE.
      TOOFEW = .FALSE.
 
*     Quick return if possible
      M = 0
.EQ.      IF( N0 ) RETURN
 
*     Simplification:
.EQ..AND..EQ..AND..EQ.      IF( IRANGE3  IL1  IUN ) IRANGE = 1
 
*     Get machine constants
      EPS = SLAMCH( 'p' )
      UFLOW = SLAMCH( 'u' )
 
 
*     Special Case when N=1
*     Treat case of 1x1 matrix for quick return
.EQ.      IF( N1 ) THEN
.EQ..OR.         IF( (IRANGE1)
.EQ..AND..GT..AND..LE..OR.     $       ((IRANGE2)(D(1)VL)(D(1)VU))
.EQ..AND..EQ..AND..EQ.     $       ((IRANGE3)(IL1)(IU1)) ) 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 )
.LE.      IF( NB1 ) 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*N - FUDGE*TWO*PIVMIN
      GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
      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
      ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN
      
.EQ.      IF( IRANGE3 ) 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 )
.NE.         IF( IINFO  0 ) THEN
            INFO = IINFO
            RETURN
         END IF
*        On exit, output intervals may not be ordered by ascending negcount
.EQ.         IF( IWORK( 6 )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.
.LT..OR..GE..OR..LT..OR..GT.         IF( NWL0  NWLN  NWU1  NWUN ) THEN
            INFO = 4
            RETURN
         END IF
 
.EQ.      ELSEIF( IRANGE2 ) THEN
         WL = VL
         WU = VU
         
.EQ.      ELSEIF( IRANGE1 ) 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
*
.EQ.         IF( IRANGE1 ) THEN
.LT..OR..GT.            IF( (IENDDOL)(IBEGINDOU) ) THEN
*              the local block contains none of eigenvalues that matter
*              to this processor
               NWU = NWU + IN
               DO 30 J = 1, IN
                  M = M + 1
                  IBLOCK( M ) = JBLK
 30            CONTINUE
               GO TO 70
            END IF
         END IF
 
.EQ.         IF( IN1 ) THEN
*           1x1 block
.GE.            IF( WLD( IBEGIN )-PIVMIN )
     $         NWL = NWL + 1
.GE.            IF( WUD( IBEGIN )-PIVMIN )
     $         NWU = NWU + 1
.EQ..OR..LT..AND..GE.            IF( IRANGE1  ( WLD( IBEGIN )-PIVMIN  WU
     $          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
         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
            SPDIAM = GU - GL            
            GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN
            GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN
*
.GT.            IF( IRANGE1 ) THEN
.LT.               IF( GUWL ) 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 )
.GE.               IF( GLGU )
     $            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 )
.NE.            IF( IINFO  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 )
.NE.            IF( IINFO  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 ) )
.GT.               IF( JIOUT-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.
.EQ.      IF( IRANGE3 ) THEN
         IDISCL = IL - 1 - NWL
         IDISCU = NWU - IU
*
.GT.         IF( IDISCL0 ) 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.
.LE..AND..GT.               IF( W( JE )WLU  IDISCL0 ) 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
.GT.         IF( IDISCU0 ) 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
.GE..AND..GT.               IF( W( JE )WUL  IDISCU0 ) 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
         
.GT..OR..GT.         IF( IDISCL0  IDISCU0 ) 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
.GT.            IF( IDISCL0 ) THEN
               WKILL = WU
               DO 100 JDISC = 1, IDISCL
                  IW = 0
                  DO 90 JE = 1, M
.NE..AND.                     IF( IBLOCK( JE )0 
.LT..OR..EQ.     $                    ( W( JE )WKILL  IW0 ) ) THEN
                        IW = JE
                        WKILL = W( JE )
                     END IF
 90               CONTINUE
                  IBLOCK( IW ) = 0
 100           CONTINUE
            END IF
.GT.            IF( IDISCU0 ) THEN
               WKILL = WL
               DO 120 JDISC = 1, IDISCU
                  IW = 0
                  DO 110 JE = 1, M
.NE..AND.                     IF( IBLOCK( JE )0 
.GE..OR..EQ.     $                    ( W( JE )WKILL  IW0 ) ) 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
.NE.               IF( IBLOCK( JE )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
.LT..OR..LT.         IF( IDISCL0  IDISCU0 ) THEN
            TOOFEW = .TRUE.
         END IF
      END IF
*
.EQ..AND..NE..OR.      IF(( IRANGE1  MN )
.EQ..AND..NE.     $   ( IRANGE3  MIU-IL+1 ) ) THEN
         TOOFEW = .TRUE.
      END IF
 
*     If ORDER='B',(IBLOCK = 2), do nothing  the eigenvalues are already sorted
*        by block.
*     If ORDER='E',(IBLOCK = 1), sort the eigenvalues from smallest to largest
      
.EQ..AND..GT.      IF( IORDER1  NSPLIT1 ) THEN
         DO 150 JE = 1, M - 1
            IE = 0
            TMP1 = W( JE )
            DO 140 J = JE + 1, M
.LT.               IF( W( J )TMP1 ) THEN
                  IE = J
                  TMP1 = W( J )
               END IF
  140       CONTINUE
.NE.            IF( IE0 ) 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 SLARRD2
*
      END