dstebz.f

*> \brief \b DSTEBZ
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DSTEBZ + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstebz.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstebz.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstebz.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
*                          M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
*                          INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          ORDER, RANGE
*       INTEGER            IL, INFO, IU, M, N, NSPLIT
*       DOUBLE PRECISION   ABSTOL, VL, VU
*       ..
*       .. Array Arguments ..
*       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
*       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DSTEBZ computes the eigenvalues of a symmetric tridiagonal
*> matrix T.  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 "Accurate Eigenvalues of a Symmetric Tridiagonal
*> Matrix", Report CS41, Computer Science Dept., Stanford
*> University, July 21, 1966.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] RANGE
*> \verbatim
*>          RANGE is CHARACTER*1
*>          = '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.
*> \endverbatim
*>
*> \param[in] ORDER
*> \verbatim
*>          ORDER is CHARACTER*1
*>          = '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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the tridiagonal matrix T.  N >= 0.
*> \endverbatim
*>
*> \param[in] VL
*> \verbatim
*>          VL is DOUBLE PRECISION
*>
*>          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'.
*> \endverbatim
*>
*> \param[in] VU
*> \verbatim
*>          VU is DOUBLE PRECISION
*>
*>          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'.
*> \endverbatim
*>
*> \param[in] IL
*> \verbatim
*>          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'.
*> \endverbatim
*>
*> \param[in] IU
*> \verbatim
*>          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'.
*> \endverbatim
*>
*> \param[in] ABSTOL
*> \verbatim
*>          ABSTOL is DOUBLE PRECISION
*>          The absolute tolerance for the eigenvalues.  An eigenvalue
*>          (or cluster) is considered to be located if it has been
*>          determined to lie in an interval whose width is ABSTOL or
*>          less.  If ABSTOL is less than or equal to zero, then ULP*|T|
*>          will be used, where |T| means the 1-norm of T.
*>
*>          Eigenvalues will be computed most accurately when ABSTOL is
*>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>          The n diagonal elements of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (N-1)
*>          The (n-1) off-diagonal elements of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[out] M
*> \verbatim
*>          M is INTEGER
*>          The actual number of eigenvalues found. 0 <= M <= N.
*>          (See also the description of INFO=2,3.)
*> \endverbatim
*>
*> \param[out] NSPLIT
*> \verbatim
*>          NSPLIT is INTEGER
*>          The number of diagonal blocks in the matrix T.
*>          1 <= NSPLIT <= N.
*> \endverbatim
*>
*> \param[out] W
*> \verbatim
*>          W is DOUBLE PRECISION array, dimension (N)
*>          On exit, the first M elements of W will contain the
*>          eigenvalues.  (DSTEBZ may use the remaining N-M elements as
*>          workspace.)
*> \endverbatim
*>
*> \param[out] IBLOCK
*> \verbatim
*>          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.  (DSTEBZ may use the remaining N-M elements as
*>          workspace.)
*> \endverbatim
*>
*> \param[out] ISPLIT
*> \verbatim
*>          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.)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (4*N)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (3*N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          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 "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.
*> \endverbatim
*
*> \par Internal Parameters:
*  =========================
*>
*> \verbatim
*>  RELFAC  DOUBLE PRECISION, default = 2.0e0
*>          The relative tolerance.  An interval (a,b] lies within
*>          "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|),
*>          where "ulp" is the machine precision (distance from 1 to
*>          the next larger floating point number.)
*>
*>  FUDGE   DOUBLE PRECISION, default = 2
*>          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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup stebz
*
*  =====================================================================
      SUBROUTINE dstebz( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
     $                   M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
     $                   INFO )
*
*  -- LAPACK computational 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
      DOUBLE PRECISION   ABSTOL, VL, VU
*     ..
*     .. Array Arguments ..
      INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO, HALF
      parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
     $                   half = 1.0d0 / two )
      DOUBLE PRECISION   FUDGE, RELFAC
      parameter( fudge = 2.1d0, relfac = 2.0d0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NCNVRG, TOOFEW
      INTEGER            IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
     $                   im, in, ioff, iorder, iout, irange, itmax,
     $                   itmp1, iw, iwoff, j, jb, jdisc, je, nb, nwl,
     $                   nwu
      DOUBLE PRECISION   ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
     $                   tmp1, tmp2, tnorm, ulp, wkill, wl, wlu, wu, wul
*     ..
*     .. Local Arrays ..
      INTEGER            IDUMMA( 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           lsame, ilaenv, dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlaebz, xerbla
*     ..
*     .. 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
*
      IF( irange.LE.0 ) THEN
         info = -1
      ELSE IF( iorder.LE.0 ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( irange.EQ.2 ) THEN
         IF( vl.GE.vu )
     $      info = -5
      ELSE IF( irange.EQ.3 .AND. ( il.LT.1 .OR. il.GT.max( 1, n ) ) )
     $          THEN
         info = -6
      ELSE IF( irange.EQ.3 .AND. ( iu.LT.min( n, il ) .OR. iu.GT.n ) )
     $          THEN
         info = -7
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DSTEBZ', -info )
         RETURN
      END IF
*
*     Initialize error flags
*
      info = 0
      ncnvrg = .false.
      toofew = .false.
*
*     Quick return if possible
*
      m = 0
      IF( n.EQ.0 )
     $   RETURN
*
*     Simplifications:
*
      IF( irange.EQ.3 .AND. il.EQ.1 .AND. iu.EQ.n )
     $   irange = 1
*
*     Get machine constants
*     NB is the minimum vector length for vector bisection, or 0
*     if only scalar is to be done.
*
      safemn = dlamch( 'S' )
      ulp = dlamch( 'P' )
      rtoli = ulp*relfac
      nb = ilaenv( 1, 'DSTEBZ', ' ', n, -1, -1, -1 )
      IF( nb.LE.1 )
     $   nb = 0
*
*     Special Case when N=1
*
      IF( n.EQ.1 ) THEN
         nsplit = 1
         isplit( 1 ) = 1
         IF( irange.EQ.2 .AND. ( vl.GE.d( 1 ) .OR. vu.LT.d( 1 ) ) ) THEN
            m = 0
         ELSE
            w( 1 ) = d( 1 )
            iblock( 1 ) = 1
            m = 1
         END IF
         RETURN
      END IF
*
*     Compute Splitting Points
*
      nsplit = 1
      work( n ) = zero
      pivmin = one
*
      DO 10 j = 2, n
         tmp1 = e( j-1 )**2
         IF( abs( d( j )*d( j-1 ) )*ulp**2+safemn.GT.tmp1 ) THEN
            isplit( nsplit ) = j - 1
            nsplit = nsplit + 1
            work( j-1 ) = zero
         ELSE
            work( j-1 ) = tmp1
            pivmin = max( pivmin, tmp1 )
         END IF
   10 CONTINUE
      isplit( nsplit ) = n
      pivmin = pivmin*safemn
*
*     Compute Interval and ATOLI
*
      IF( irange.EQ.3 ) THEN
*
*        RANGE='I': Compute the interval containing eigenvalues
*                   IL through IU.
*
*        Compute Gershgorin interval for entire (split) matrix
*        and use it as the initial interval
*
         gu = d( 1 )
         gl = d( 1 )
         tmp1 = zero
*
         DO 20 j = 1, n - 1
            tmp2 = sqrt( work( j ) )
            gu = max( gu, d( j )+tmp1+tmp2 )
            gl = min( gl, d( j )-tmp1-tmp2 )
            tmp1 = tmp2
   20    CONTINUE
*
         gu = max( gu, d( n )+tmp1 )
         gl = min( gl, d( n )-tmp1 )
         tnorm = max( abs( gl ), abs( gu ) )
         gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
         gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
*
*        Compute Iteration parameters
*
         itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
     $           log( two ) ) + 2
         IF( abstol.LE.zero ) THEN
            atoli = ulp*tnorm
         ELSE
            atoli = abstol
         END IF
*
         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 dlaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin, d, e,
     $                work, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
     $                iwork, w, iblock, iinfo )
*
         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
*
         IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
            info = 4
            RETURN
         END IF
      ELSE
*
*        RANGE='A' or 'V' -- Set ATOLI
*
         tnorm = max( abs( d( 1 ) )+abs( e( 1 ) ),
     $           abs( d( n ) )+abs( e( n-1 ) ) )
*
         DO 30 j = 2, n - 1
            tnorm = max( tnorm, abs( d( j ) )+abs( e( j-1 ) )+
     $              abs( e( j ) ) )
   30    CONTINUE
*
         IF( abstol.LE.zero ) THEN
            atoli = ulp*tnorm
         ELSE
            atoli = abstol
         END IF
*
         IF( irange.EQ.2 ) THEN
            wl = vl
            wu = vu
         ELSE
            wl = zero
            wu = zero
         END IF
      END IF
*
*     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 jb = 1, nsplit
         ioff = iend
         ibegin = ioff + 1
         iend = isplit( jb )
         in = iend - ioff
*
         IF( in.EQ.1 ) THEN
*
*           Special Case -- IN=1
*
            IF( irange.EQ.1 .OR. wl.GE.d( ibegin )-pivmin )
     $         nwl = nwl + 1
            IF( irange.EQ.1 .OR. wu.GE.d( ibegin )-pivmin )
     $         nwu = nwu + 1
            IF( irange.EQ.1 .OR. ( wl.LT.d( ibegin )-pivmin .AND. wu.GE.
     $          d( ibegin )-pivmin ) ) THEN
               m = m + 1
               w( m ) = d( ibegin )
               iblock( m ) = jb
            END IF
         ELSE
*
*           General Case -- IN > 1
*
*           Compute Gershgorin Interval
*           and use it as the initial interval
*
            gu = d( ibegin )
            gl = d( ibegin )
            tmp1 = zero
*
            DO 40 j = ibegin, iend - 1
               tmp2 = abs( e( j ) )
               gu = max( gu, d( j )+tmp1+tmp2 )
               gl = min( gl, d( j )-tmp1-tmp2 )
               tmp1 = tmp2
   40       CONTINUE
*
            gu = max( gu, d( iend )+tmp1 )
            gl = min( gl, d( iend )-tmp1 )
            bnorm = max( abs( gl ), abs( gu ) )
            gl = gl - fudge*bnorm*ulp*in - fudge*pivmin
            gu = gu + fudge*bnorm*ulp*in + fudge*pivmin
*
*           Compute ATOLI for the current submatrix
*
            IF( abstol.LE.zero ) THEN
               atoli = ulp*max( abs( gl ), abs( gu ) )
            ELSE
               atoli = abstol
            END IF
*
            IF( irange.GT.1 ) THEN
               IF( gu.LT.wl ) THEN
                  nwl = nwl + in
                  nwu = nwu + in
                  GO TO 70
               END IF
               gl = max( gl, wl )
               gu = min( gu, wu )
               IF( gl.GE.gu )
     $            GO TO 70
            END IF
*
*           Set Up Initial Interval
*
            work( n+1 ) = gl
            work( n+in+1 ) = gu
            CALL dlaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
     $                   d( ibegin ), e( ibegin ), work( ibegin ),
     $                   idumma, work( n+1 ), work( n+2*in+1 ), im,
     $                   iwork, w( m+1 ), iblock( m+1 ), iinfo )
*
            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 dlaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
     $                   d( ibegin ), e( ibegin ), work( ibegin ),
     $                   idumma, work( n+1 ), work( n+2*in+1 ), iout,
     $                   iwork, w( m+1 ), iblock( m+1 ), iinfo )
*
*           Copy Eigenvalues Into W and IBLOCK
*           Use -JB for block number for unconverged eigenvalues.
*
            DO 60 j = 1, iout
               tmp1 = half*( work( j+n )+work( j+in+n ) )
*
*              Flag non-convergence.
*
               IF( j.GT.iout-iinfo ) THEN
                  ncnvrg = .true.
                  ib = -jb
               ELSE
                  ib = jb
               END IF
               DO 50 je = iwork( j ) + 1 + iwoff,
     $                 iwork( j+in ) + iwoff
                  w( je ) = tmp1
                  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.3 ) THEN
         im = 0
         idiscl = il - 1 - nwl
         idiscu = nwu - iu
*
         IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
            DO 80 je = 1, m
               IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
                  idiscl = idiscl - 1
               ELSE IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
                  idiscu = idiscu - 1
               ELSE
                  im = im + 1
                  w( im ) = w( je )
                  iblock( im ) = iblock( je )
               END IF
   80       CONTINUE
            m = im
         END IF
         IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
*
*           Code to deal with effects of bad arithmetic:
*           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 simply finding the smallest/largest
*           eigenvalue(s).
*
*           (If N(w) is monotone non-decreasing, this should never
*               happen.)
*
            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 ).GT.wkill .OR. iw.EQ.0 ) ) THEN
                        iw = je
                        wkill = w( je )
                     END IF
  110             CONTINUE
                  iblock( iw ) = 0
  120          CONTINUE
            END IF
            im = 0
            DO 130 je = 1, m
               IF( iblock( je ).NE.0 ) THEN
                  im = im + 1
                  w( im ) = w( 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 ORDER='B', do nothing -- the eigenvalues are already sorted
*        by block.
*     If ORDER='E', sort the eigenvalues from smallest to largest
*
      IF( iorder.EQ.1 .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
               itmp1 = iblock( ie )
               w( ie ) = w( je )
               iblock( ie ) = iblock( je )
               w( je ) = tmp1
               iblock( je ) = itmp1
            END IF
  150    CONTINUE
      END IF
*
      info = 0
      IF( ncnvrg )
     $   info = info + 1
      IF( toofew )
     $   info = info + 2
      RETURN
*
*     End of DSTEBZ
*
      END