```      SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
\$                   M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
\$                   INFO )
*
*  -- LAPACK routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*     8-18-00:  Increase FUDGE factor for T3E (eca)
*
*     .. 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( * )
*     ..
*
*  Purpose
*  =======
*
*  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.
*
*  Arguments
*  =========
*
*  RANGE   (input) 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.
*
*  ORDER   (input) 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.
*
*  N       (input) INTEGER
*          The order of the tridiagonal matrix T.  N >= 0.
*
*  VL      (input) DOUBLE PRECISION
*  VU      (input) DOUBLE PRECISION
*          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'.
*
*  ABSTOL  (input) 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.
*
*  D       (input) DOUBLE PRECISION array, dimension (N)
*          The n diagonal elements of the tridiagonal matrix T.
*
*  E       (input) DOUBLE PRECISION array, dimension (N-1)
*          The (n-1) off-diagonal elements of the tridiagonal matrix T.
*
*  M       (output) INTEGER
*          The actual number of eigenvalues found. 0 <= M <= N.
*
*  NSPLIT  (output) INTEGER
*          The number of diagonal blocks in the matrix T.
*          1 <= NSPLIT <= N.
*
*  W       (output) 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.)
*
*  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.  (DSTEBZ may use the remaining N-M elements as
*          workspace.)
*
*  ISPLIT  (output) 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.)
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
*
*  IWORK   (workspace) INTEGER array, dimension (3*N)
*
*  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
*  ===================
*
*  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.
*
*  =====================================================================
*
*     .. 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
*
*DIR\$ NOVECTOR
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

```