SUBROUTINE PDSTEBZ( ICTXT, RANGE, ORDER, N, VL, VU, IL, IU,
     $                    ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT,
     $                    WORK, LWORK, IWORK, LIWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*     .. Scalar Arguments ..
      CHARACTER          ORDER, RANGE
      INTEGER            ICTXT, IL, INFO, IU, LIWORK, LWORK, M, N,
     $                   NSPLIT
      DOUBLE PRECISION   ABSTOL, VL, VU
*     ..
*     .. Array Arguments ..
      INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PDSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix in
*  parallel. The user may ask for all eigenvalues, all eigenvalues in
*  the interval [VL, VU], or the eigenvalues indexed IL through IU. A
*  static partitioning of work is done at the beginning of PDSTEBZ which
*  results in all processes finding an (almost) equal number of
*  eigenvalues.
*
*  NOTE : It is assumed that the user is on an IEEE machine. If the user
*         is not on an IEEE mchine, set the compile time flag NO_IEEE
*         to 1 (in SLmake.inc). The features of IEEE arithmetic that
*         are needed for the "fast" Sturm Count are : (a) infinity
*         arithmetic (b) the sign bit of a single precision floating
*         point number is assumed be in the 32nd bit position
*         (c) the sign of negative zero.
*
*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
*  Matrix", Report CS41, Computer Science Dept., Stanford
*  University, July 21, 1966.
*
*  Arguments
*  =========
*
*  ICTXT   (global input) INTEGER
*          The BLACS context handle.
*
*  RANGE   (global input) CHARACTER
*          Specifies which eigenvalues are to be found.
*          = 'A': ("All")   all eigenvalues will be found.
*          = 'V': ("Value") all eigenvalues in the interval
*                           [VL, VU] will be found.
*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
*                           entire matrix) will be found.
*
*  ORDER   (global input) CHARACTER
*          Specifies the order in which the eigenvalues and their block
*          numbers are stored in W and IBLOCK.
*          = '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       (global input) INTEGER
*          The order of the tridiagonal matrix T.  N >= 0.
*
*  VL      (global input) DOUBLE PRECISION
*          If RANGE='V', the lower bound of the interval to be searched
*          for eigenvalues.  Eigenvalues less than VL will not be
*          returned.  Not referenced if RANGE='A' or 'I'.
*
*  VU      (global input) DOUBLE PRECISION
*          If RANGE='V', the upper bound of the interval to be searched
*          for eigenvalues.  Eigenvalues greater than VU will not be
*          returned.  VU must be greater than VL.  Not referenced if
*          RANGE='A' or 'I'.
*
*  IL      (global input) INTEGER
*          If RANGE='I', the index (from smallest to largest) of the
*          smallest eigenvalue to be returned.  IL must be at least 1.
*          Not referenced if RANGE='A' or 'V'.
*
*  IU      (global input) INTEGER
*          If RANGE='I', the index (from smallest to largest) of the
*          largest eigenvalue to be returned.  IU must be at least IL
*          and no greater than N.  Not referenced if RANGE='A' or 'V'.
*
*  ABSTOL  (global 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 the underflow threshold DLAMCH('U'), not zero.
*          Note : If eigenvectors are desired later by inverse iteration
*          ( PDSTEIN ), ABSTOL should be set to 2*PDLAMCH('S').
*
*  D       (global input) DOUBLE PRECISION array, dimension (N)
*          The n diagonal elements of the tridiagonal matrix T.  To
*          avoid overflow, the matrix must be scaled so that its largest
*          entry 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.
*
*  E       (global input) DOUBLE PRECISION array, dimension (N-1)
*          The (n-1) off-diagonal elements of the tridiagonal matrix T.
*          To avoid overflow, the matrix must be scaled so that its
*          largest entry 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.
*
*  M       (global output) INTEGER
*          The actual number of eigenvalues found. 0 <= M <= N.
*          (See also the description of INFO=2)
*
*  NSPLIT  (global output) INTEGER
*          The number of diagonal blocks in the matrix T.
*          1 <= NSPLIT <= N.
*
*  W       (global output) DOUBLE PRECISION array, dimension (N)
*          On exit, the first M elements of W contain the eigenvalues
*          on all processes.
*
*  IBLOCK  (global 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 IBLOCK(i) specifies which block (from 1
*          to the number of blocks) the eigenvalue W(i) belongs to.
*          NOTE:  in the (theoretically impossible) event that bisection
*          does not converge for some or all eigenvalues, INFO is set
*          to 1 and the ones for which it did not are identified by a
*          negative block number.
*
*  ISPLIT  (global 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    (local workspace) DOUBLE PRECISION array,
*          dimension ( MAX( 5*N, 7 ) )
*
*  LWORK   (local input) INTEGER
*          size of array WORK must be >= MAX( 5*N, 7 )
*          If LWORK = -1, then LWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  IWORK   (local workspace) INTEGER array, dimension ( MAX( 4*N, 14 ) )
*
*  LIWORK  (local input) INTEGER
*          size of array IWORK must be >= MAX( 4*N, 14, NPROCS )
*          If LIWORK = -1, then LIWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  INFO    (global 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 : 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 arithmetic
*                    which is less accurate than PDLAMCH says.
*              = 2 : There is a mismatch between the number of
*                    eigenvalues output and the number desired.
*              = 3 : RANGE='i', and the Gershgorin interval initially
*                    used was incorrect. 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.0
*          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.0
*          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 the accuracy of the solution.
*
*  =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, ICHAR, MAX, MIN, MOD
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            BLACS_PNUM
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           LSAME, BLACS_PNUM, PDLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_FREEBUFF, BLACS_GET, BLACS_GRIDEXIT,
     $                   BLACS_GRIDINFO, BLACS_GRIDMAP, DGEBR2D,
     $                   DGEBS2D, DGERV2D, DGESD2D, DLASRT2, GLOBCHK,
     $                   IGEBR2D, IGEBS2D, IGERV2D, IGESD2D, IGSUM2D,
     $                   PDLAEBZ, PDLAIECTB, PDLAIECTL, PDLAPDCT,
     $                   PDLASNBT, PXERBLA
*     ..
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
     $                   MB_, NB_, RSRC_, CSRC_, LLD_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      INTEGER            BIGNUM, DESCMULT
      PARAMETER          ( BIGNUM = 10000, DESCMULT = 100 )
      DOUBLE PRECISION   ZERO, ONE, TWO, FIVE, HALF
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
     $                   FIVE = 5.0D+0, HALF = 1.0D+0 / TWO )
      DOUBLE PRECISION   FUDGE, RELFAC
      PARAMETER          ( FUDGE = 2.0D+0, RELFAC = 2.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY
      INTEGER            BLKNO, FOUND, I, IBEGIN, IEFLAG, IEND, IFRST,
     $                   IINFO, ILAST, ILOAD, IM, IMYLOAD, IN, INDRIW1,
     $                   INDRIW2, INDRW1, INDRW2, INXTLOAD, IOFF,
     $                   IORDER, IOUT, IRANGE, IRECV, IREM, ITMP1,
     $                   ITMP2, J, JB, K, LAST, LEXTRA, LREQ, MYCOL,
     $                   MYROW, NALPHA, NBETA, NCMP, NEIGINT, NEXT, NGL,
     $                   NGLOB, NGU, NINT, NPCOL, NPROW, OFFSET,
     $                   ONEDCONTEXT, P, PREV, REXTRA, RREQ, SELF,
     $                   TORECV
      DOUBLE PRECISION   ALPHA, ATOLI, BETA, BNORM, DRECV, DSEND, GL,
     $                   GU, INITVL, INITVU, LSAVE, MID, PIVMIN, RELTOL,
     $                   SAFEMN, TMP1, TMP2, TNORM, ULP
*     ..
*     .. Local Arrays ..
      INTEGER            IDUM( 5, 2 )
*     ..
*     .. Executable Statements ..
*       This is just to keep ftnchek happy
      IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
     $    RSRC_.LT.0 )RETURN
*
*     Set up process grid
*
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
      INFO = 0
      M = 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' ) .OR. LSAME( ORDER, 'A' ) ) THEN
         IORDER = 1
      ELSE
         IORDER = 0
      END IF
*
*     Check for Errors
*
      IF( NPROW.EQ.-1 ) THEN
         INFO = -1
      ELSE
*
*     Get machine constants
*
         SAFEMN = PDLAMCH( ICTXT, 'S' )
         ULP = PDLAMCH( ICTXT, 'P' )
         RELTOL = ULP*RELFAC
         IDUM( 1, 1 ) = ICHAR( RANGE )
         IDUM( 1, 2 ) = 2
         IDUM( 2, 1 ) = ICHAR( ORDER )
         IDUM( 2, 2 ) = 3
         IDUM( 3, 1 ) = N
         IDUM( 3, 2 ) = 4
         NGLOB = 5
         IF( IRANGE.EQ.3 ) THEN
            IDUM( 4, 1 ) = IL
            IDUM( 4, 2 ) = 7
            IDUM( 5, 1 ) = IU
            IDUM( 5, 2 ) = 8
         ELSE
            IDUM( 4, 1 ) = 0
            IDUM( 4, 2 ) = 0
            IDUM( 5, 1 ) = 0
            IDUM( 5, 2 ) = 0
         END IF
         IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
            WORK( 1 ) = ABSTOL
            IF( IRANGE.EQ.2 ) THEN
               WORK( 2 ) = VL
               WORK( 3 ) = VU
            ELSE
               WORK( 2 ) = ZERO
               WORK( 3 ) = ZERO
            END IF
            CALL DGEBS2D( ICTXT, 'ALL', ' ', 3, 1, WORK, 3 )
         ELSE
            CALL DGEBR2D( ICTXT, 'ALL', ' ', 3, 1, WORK, 3, 0, 0 )
         END IF
         LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
         IF( INFO.EQ.0 ) THEN
            IF( IRANGE.EQ.0 ) THEN
               INFO = -2
            ELSE IF( IORDER.EQ.0 ) THEN
               INFO = -3
            ELSE IF( IRANGE.EQ.2 .AND. VL.GE.VU ) THEN
               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
            ELSE IF( LWORK.LT.MAX( 5*N, 7 ) .AND. .NOT.LQUERY ) THEN
               INFO = -18
            ELSE IF( LIWORK.LT.MAX( 4*N, 14, NPROW*NPCOL ) .AND. .NOT.
     $              LQUERY ) THEN
               INFO = -20
            ELSE IF( IRANGE.EQ.2 .AND. ( ABS( WORK( 2 )-VL ).GT.FIVE*
     $              ULP*ABS( VL ) ) ) THEN
               INFO = -5
            ELSE IF( IRANGE.EQ.2 .AND. ( ABS( WORK( 3 )-VU ).GT.FIVE*
     $              ULP*ABS( VU ) ) ) THEN
               INFO = -6
            ELSE IF( ABS( WORK( 1 )-ABSTOL ).GT.FIVE*ULP*ABS( ABSTOL ) )
     $               THEN
               INFO = -9
            END IF
         END IF
         IF( INFO.EQ.0 )
     $      INFO = BIGNUM
         CALL GLOBCHK( ICTXT, NGLOB, IDUM, 5, IWORK, INFO )
         IF( INFO.EQ.BIGNUM ) THEN
            INFO = 0
         ELSE IF( MOD( INFO, DESCMULT ).EQ.0 ) THEN
            INFO = -INFO / DESCMULT
         ELSE
            INFO = -INFO
         END IF
      END IF
      WORK( 1 ) = DBLE( MAX( 5*N, 7 ) )
      IWORK( 1 ) = MAX( 4*N, 14, NPROW*NPCOL )
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PDSTEBZ', -INFO )
         RETURN
      ELSE IF( LWORK.EQ.-1 .AND. LIWORK.EQ.-1 ) THEN
         RETURN
      END IF
*
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
*
      K = 1
      DO 20 I = 0, NPROW - 1
         DO 10 J = 0, NPCOL - 1
            IWORK( K ) = BLACS_PNUM( ICTXT, I, J )
            K = K + 1
   10    CONTINUE
   20 CONTINUE
*
      P = NPROW*NPCOL
      NPROW = 1
      NPCOL = P
*
      CALL BLACS_GET( ICTXT, 10, ONEDCONTEXT )
      CALL BLACS_GRIDMAP( ONEDCONTEXT, IWORK, NPROW, NPROW, NPCOL )
      CALL BLACS_GRIDINFO( ONEDCONTEXT, I, J, K, SELF )
*
*     Simplifications:
*
      IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
     $   IRANGE = 1
*
      NEXT = MOD( SELF+1, P )
      PREV = MOD( P+SELF-1, P )
*
*     Compute squares of off-diagonals, splitting points and pivmin.
*     Interleave diagonals and off-diagonals.
*
      INDRW1 = MAX( 2*N, 4 )
      INDRW2 = INDRW1 + 2*N
      INDRIW1 = MAX( 2*N, 8 )
      NSPLIT = 1
      WORK( INDRW1+2*N ) = ZERO
      PIVMIN = ONE
*
      DO 30 I = 1, N - 1
         TMP1 = E( I )**2
         J = 2*I
         WORK( INDRW1+J-1 ) = D( I )
         IF( ABS( D( I+1 )*D( I ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
            ISPLIT( NSPLIT ) = I
            NSPLIT = NSPLIT + 1
            WORK( INDRW1+J ) = ZERO
         ELSE
            WORK( INDRW1+J ) = TMP1
            PIVMIN = MAX( PIVMIN, TMP1 )
         END IF
   30 CONTINUE
      WORK( INDRW1+2*N-1 ) = D( N )
      ISPLIT( NSPLIT ) = N
      PIVMIN = PIVMIN*SAFEMN
*
*     Compute Gershgorin interval [gl,gu] for entire matrix
*
      GU = D( 1 )
      GL = D( 1 )
      TMP1 = ZERO
*
      DO 40 I = 1, N - 1
         TMP2 = ABS( E( I ) )
         GU = MAX( GU, D( I )+TMP1+TMP2 )
         GL = MIN( GL, D( I )-TMP1-TMP2 )
         TMP1 = TMP2
   40 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
*
      IF( ABSTOL.LE.ZERO ) THEN
         ATOLI = ULP*TNORM
      ELSE
         ATOLI = ABSTOL
      END IF
*
*     Find out if on an IEEE machine, the sign bit is the
*     32nd bit (Big Endian) or the 64th bit (Little Endian)
*
      IF( IRANGE.EQ.1 .OR. NSPLIT.EQ.1 ) THEN
         CALL PDLASNBT( IEFLAG )
      ELSE
         IEFLAG = 0
      END IF
      LEXTRA = 0
      REXTRA = 0
*
*     Form Initial Interval containing desired eigenvalues
*
      IF( IRANGE.EQ.1 ) THEN
         INITVL = GL
         INITVU = GU
         WORK( 1 ) = GL
         WORK( 2 ) = GU
         IWORK( 1 ) = 0
         IWORK( 2 ) = N
         IFRST = 1
         ILAST = N
      ELSE IF( IRANGE.EQ.2 ) THEN
         IF( VL.GT.GL ) THEN
            IF( IEFLAG.EQ.0 ) THEN
               CALL PDLAPDCT( VL, N, WORK( INDRW1+1 ), PIVMIN, IFRST )
            ELSE IF( IEFLAG.EQ.1 ) THEN
               CALL PDLAIECTB( VL, N, WORK( INDRW1+1 ), IFRST )
            ELSE
               CALL PDLAIECTL( VL, N, WORK( INDRW1+1 ), IFRST )
            END IF
            IFRST = IFRST + 1
            INITVL = VL
         ELSE
            INITVL = GL
            IFRST = 1
         END IF
         IF( VU.LT.GU ) THEN
            IF( IEFLAG.EQ.0 ) THEN
               CALL PDLAPDCT( VU, N, WORK( INDRW1+1 ), PIVMIN, ILAST )
            ELSE IF( IEFLAG.EQ.1 ) THEN
               CALL PDLAIECTB( VU, N, WORK( INDRW1+1 ), ILAST )
            ELSE
               CALL PDLAIECTL( VU, N, WORK( INDRW1+1 ), ILAST )
            END IF
            INITVU = VU
         ELSE
            INITVU = GU
            ILAST = N
         END IF
         WORK( 1 ) = INITVL
         WORK( 2 ) = INITVU
         IWORK( 1 ) = IFRST - 1
         IWORK( 2 ) = ILAST
      ELSE IF( IRANGE.EQ.3 ) THEN
         WORK( 1 ) = GL
         WORK( 2 ) = GU
         IWORK( 1 ) = 0
         IWORK( 2 ) = N
         IWORK( 5 ) = IL - 1
         IWORK( 6 ) = IU
         CALL PDLAEBZ( 0, N, 2, 1, ATOLI, RELTOL, PIVMIN,
     $                 WORK( INDRW1+1 ), IWORK( 5 ), WORK, IWORK, NINT,
     $                 LSAVE, IEFLAG, IINFO )
         IF( IINFO.NE.0 ) THEN
            INFO = 3
            GO TO 230
         END IF
         IF( NINT.GT.1 ) THEN
            IF( IWORK( 5 ).EQ.IL-1 ) THEN
               WORK( 2 ) = WORK( 4 )
               IWORK( 2 ) = IWORK( 4 )
            ELSE
               WORK( 1 ) = WORK( 3 )
               IWORK( 1 ) = IWORK( 3 )
            END IF
            IF( IWORK( 1 ).LT.0 .OR. IWORK( 1 ).GT.IL-1 .OR.
     $          IWORK( 2 ).LE.MIN( IU-1, IWORK( 1 ) ) .OR.
     $          IWORK( 2 ).GT.N ) THEN
               INFO = 3
               GO TO 230
            END IF
         END IF
         LEXTRA = IL - 1 - IWORK( 1 )
         REXTRA = IWORK( 2 ) - IU
         INITVL = WORK( 1 )
         INITVU = WORK( 2 )
         IFRST = IL
         ILAST = IU
      END IF
*     NVL = IFRST - 1
*     NVU = ILAST
      GL = INITVL
      GU = INITVU
      NGL = IWORK( 1 )
      NGU = IWORK( 2 )
      IM = 0
      FOUND = 0
      INDRIW2 = INDRIW1 + NGU - NGL
      IEND = 0
      IF( IFRST.GT.ILAST )
     $   GO TO 100
      IF( IFRST.EQ.1 .AND. ILAST.EQ.N )
     $   IRANGE = 1
*
*     Find Eigenvalues -- Loop Over Blocks
*
      DO 90 JB = 1, NSPLIT
         IOFF = IEND
         IBEGIN = IOFF + 1
         IEND = ISPLIT( JB )
         IN = IEND - IOFF
         IF( JB.NE.1 ) THEN
            IF( IRANGE.NE.1 ) THEN
               FOUND = IM
*
*              Find total number of eigenvalues found thus far
*
               CALL IGSUM2D( ONEDCONTEXT, 'All', ' ', 1, 1, FOUND, 1,
     $                       -1, -1 )
            ELSE
               FOUND = IOFF
            END IF
         END IF
*         IF( SELF.GE.P )
*     $      GO TO 30
         IF( IN.NE.N ) THEN
*
*           Compute Gershgorin interval [gl,gu] for split matrix
*
            GU = D( IBEGIN )
            GL = D( IBEGIN )
            TMP1 = ZERO
*
            DO 50 J = IBEGIN, IEND - 1
               TMP2 = ABS( E( J ) )
               GU = MAX( GU, D( J )+TMP1+TMP2 )
               GL = MIN( GL, D( J )-TMP1-TMP2 )
               TMP1 = TMP2
   50       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*BNORM
            ELSE
               ATOLI = ABSTOL
            END IF
*
            IF( GL.LT.INITVL ) THEN
               GL = INITVL
               IF( IEFLAG.EQ.0 ) THEN
                  CALL PDLAPDCT( GL, IN, WORK( INDRW1+2*IOFF+1 ),
     $                           PIVMIN, NGL )
               ELSE IF( IEFLAG.EQ.1 ) THEN
                  CALL PDLAIECTB( GL, IN, WORK( INDRW1+2*IOFF+1 ), NGL )
               ELSE
                  CALL PDLAIECTL( GL, IN, WORK( INDRW1+2*IOFF+1 ), NGL )
               END IF
            ELSE
               NGL = 0
            END IF
            IF( GU.GT.INITVU ) THEN
               GU = INITVU
               IF( IEFLAG.EQ.0 ) THEN
                  CALL PDLAPDCT( GU, IN, WORK( INDRW1+2*IOFF+1 ),
     $                           PIVMIN, NGU )
               ELSE IF( IEFLAG.EQ.1 ) THEN
                  CALL PDLAIECTB( GU, IN, WORK( INDRW1+2*IOFF+1 ), NGU )
               ELSE
                  CALL PDLAIECTL( GU, IN, WORK( INDRW1+2*IOFF+1 ), NGU )
               END IF
            ELSE
               NGU = IN
            END IF
            IF( NGL.GE.NGU )
     $         GO TO 90
            WORK( 1 ) = GL
            WORK( 2 ) = GU
            IWORK( 1 ) = NGL
            IWORK( 2 ) = NGU
         END IF
         OFFSET = FOUND - NGL
         BLKNO = JB
*
*        Do a static partitioning of work so that each process
*        has to find an (almost) equal number of eigenvalues
*
         NCMP = NGU - NGL
         ILOAD = NCMP / P
         IREM = NCMP - ILOAD*P
         ITMP1 = MOD( SELF-FOUND, P )
         IF( ITMP1.LT.0 )
     $      ITMP1 = ITMP1 + P
         IF( ITMP1.LT.IREM ) THEN
            IMYLOAD = ILOAD + 1
         ELSE
            IMYLOAD = ILOAD
         END IF
         IF( IMYLOAD.EQ.0 ) THEN
            GO TO 90
         ELSE IF( IN.EQ.1 ) THEN
            WORK( INDRW2+IM+1 ) = WORK( INDRW1+2*IOFF+1 )
            IWORK( INDRIW1+IM+1 ) = BLKNO
            IWORK( INDRIW2+IM+1 ) = OFFSET + 1
            IM = IM + 1
            GO TO 90
         ELSE
            INXTLOAD = ILOAD
            ITMP2 = MOD( SELF+1-FOUND, P )
            IF( ITMP2.LT.0 )
     $         ITMP2 = ITMP2 + P
            IF( ITMP2.LT.IREM )
     $         INXTLOAD = INXTLOAD + 1
            LREQ = NGL + ITMP1*ILOAD + MIN( IREM, ITMP1 )
            RREQ = LREQ + IMYLOAD
            IWORK( 5 ) = LREQ
            IWORK( 6 ) = RREQ
            TMP1 = WORK( 1 )
            ITMP1 = IWORK( 1 )
            CALL PDLAEBZ( 1, IN, 1, 1, ATOLI, RELTOL, PIVMIN,
     $                    WORK( INDRW1+2*IOFF+1 ), IWORK( 5 ), WORK,
     $                    IWORK, NINT, LSAVE, IEFLAG, IINFO )
            ALPHA = WORK( 1 )
            BETA = WORK( 2 )
            NALPHA = IWORK( 1 )
            NBETA = IWORK( 2 )
            DSEND = BETA
            IF( NBETA.GT.RREQ+INXTLOAD ) THEN
               NBETA = RREQ
               DSEND = ALPHA
            END IF
            LAST = MOD( FOUND+MIN( NGU-NGL, P )-1, P )
            IF( LAST.LT.0 )
     $         LAST = LAST + P
            IF( SELF.NE.LAST ) THEN
               CALL DGESD2D( ONEDCONTEXT, 1, 1, DSEND, 1, 0, NEXT )
               CALL IGESD2D( ONEDCONTEXT, 1, 1, NBETA, 1, 0, NEXT )
            END IF
            IF( SELF.NE.MOD( FOUND, P ) ) THEN
               CALL DGERV2D( ONEDCONTEXT, 1, 1, DRECV, 1, 0, PREV )
               CALL IGERV2D( ONEDCONTEXT, 1, 1, IRECV, 1, 0, PREV )
            ELSE
               DRECV = TMP1
               IRECV = ITMP1
            END IF
            WORK( 1 ) = MAX( LSAVE, DRECV )
            IWORK( 1 ) = IRECV
            ALPHA = MAX( ALPHA, WORK( 1 ) )
            NALPHA = MAX( NALPHA, IRECV )
            IF( BETA-ALPHA.LE.MAX( ATOLI, RELTOL*MAX( ABS( ALPHA ),
     $          ABS( BETA ) ) ) ) THEN
               MID = HALF*( ALPHA+BETA )
               DO 60 J = OFFSET + NALPHA + 1, OFFSET + NBETA
                  WORK( INDRW2+IM+1 ) = MID
                  IWORK( INDRIW1+IM+1 ) = BLKNO
                  IWORK( INDRIW2+IM+1 ) = J
                  IM = IM + 1
   60          CONTINUE
               WORK( 2 ) = ALPHA
               IWORK( 2 ) = NALPHA
            END IF
         END IF
         NEIGINT = IWORK( 2 ) - IWORK( 1 )
         IF( NEIGINT.LE.0 )
     $      GO TO 90
*
*        Call the main computational routine
*
         CALL PDLAEBZ( 2, IN, NEIGINT, 1, ATOLI, RELTOL, PIVMIN,
     $                 WORK( INDRW1+2*IOFF+1 ), IWORK, WORK, IWORK,
     $                 IOUT, LSAVE, IEFLAG, IINFO )
         IF( IINFO.NE.0 ) THEN
            INFO = 1
         END IF
         DO 80 I = 1, IOUT
            MID = HALF*( WORK( 2*I-1 )+WORK( 2*I ) )
            IF( I.GT.IOUT-IINFO )
     $         BLKNO = -BLKNO
            DO 70 J = OFFSET + IWORK( 2*I-1 ) + 1,
     $              OFFSET + IWORK( 2*I )
               WORK( INDRW2+IM+1 ) = MID
               IWORK( INDRIW1+IM+1 ) = BLKNO
               IWORK( INDRIW2+IM+1 ) = J
               IM = IM + 1
   70       CONTINUE
   80    CONTINUE
   90 CONTINUE
*
*     Find out total number of eigenvalues computed
*
  100 CONTINUE
      M = IM
      CALL IGSUM2D( ONEDCONTEXT, 'ALL', ' ', 1, 1, M, 1, -1, -1 )
*
*     Move the eigenvalues found to their final destinations
*
      DO 130 I = 1, P
         IF( SELF.EQ.I-1 ) THEN
            CALL IGEBS2D( ONEDCONTEXT, 'ALL', ' ', 1, 1, IM, 1 )
            IF( IM.NE.0 ) THEN
               CALL IGEBS2D( ONEDCONTEXT, 'ALL', ' ', IM, 1,
     $                       IWORK( INDRIW2+1 ), IM )
               CALL DGEBS2D( ONEDCONTEXT, 'ALL', ' ', IM, 1,
     $                       WORK( INDRW2+1 ), IM )
               CALL IGEBS2D( ONEDCONTEXT, 'ALL', ' ', IM, 1,
     $                       IWORK( INDRIW1+1 ), IM )
               DO 110 J = 1, IM
                  W( IWORK( INDRIW2+J ) ) = WORK( INDRW2+J )
                  IBLOCK( IWORK( INDRIW2+J ) ) = IWORK( INDRIW1+J )
  110          CONTINUE
            END IF
         ELSE
            CALL IGEBR2D( ONEDCONTEXT, 'ALL', ' ', 1, 1, TORECV, 1, 0,
     $                    I-1 )
            IF( TORECV.NE.0 ) THEN
               CALL IGEBR2D( ONEDCONTEXT, 'ALL', ' ', TORECV, 1, IWORK,
     $                       TORECV, 0, I-1 )
               CALL DGEBR2D( ONEDCONTEXT, 'ALL', ' ', TORECV, 1, WORK,
     $                       TORECV, 0, I-1 )
               CALL IGEBR2D( ONEDCONTEXT, 'ALL', ' ', TORECV, 1,
     $                       IWORK( N+1 ), TORECV, 0, I-1 )
               DO 120 J = 1, TORECV
                  W( IWORK( J ) ) = WORK( J )
                  IBLOCK( IWORK( J ) ) = IWORK( N+J )
  120          CONTINUE
            END IF
         END IF
  130 CONTINUE
      IF( NSPLIT.GT.1 .AND. IORDER.EQ.1 ) THEN
*
*        Sort the eigenvalues
*
*
         DO 140 I = 1, M
            IWORK( M+I ) = I
  140    CONTINUE
         CALL DLASRT2( 'I', M, W, IWORK( M+1 ), IINFO )
         DO 150 I = 1, M
            IWORK( I ) = IBLOCK( I )
  150    CONTINUE
         DO 160 I = 1, M
            IBLOCK( I ) = IWORK( IWORK( M+I ) )
  160    CONTINUE
      END IF
      IF( IRANGE.EQ.3 .AND. ( LEXTRA.GT.0 .OR. REXTRA.GT.0 ) ) THEN
*
*        Discard unwanted eigenvalues (occurs only when RANGE = 'I',
*        and eigenvalues IL, and/or IU are in a cluster)
*
         DO 170 I = 1, M
            WORK( I ) = W( I )
            IWORK( I ) = I
            IWORK( M+I ) = I
  170    CONTINUE
         DO 190 I = 1, LEXTRA
            ITMP1 = I
            DO 180 J = I + 1, M
               IF( WORK( J ).LT.WORK( ITMP1 ) ) THEN
                  ITMP1 = J
               END IF
  180       CONTINUE
            TMP1 = WORK( I )
            WORK( I ) = WORK( ITMP1 )
            WORK( ITMP1 ) = TMP1
            IWORK( IWORK( M+ITMP1 ) ) = I
            IWORK( IWORK( M+I ) ) = ITMP1
            ITMP2 = IWORK( M+I )
            IWORK( M+I ) = IWORK( M+ITMP1 )
            IWORK( M+ITMP1 ) = ITMP2
  190    CONTINUE
         DO 210 I = 1, REXTRA
            ITMP1 = M - I + 1
            DO 200 J = M - I, LEXTRA + 1, -1
               IF( WORK( J ).GT.WORK( ITMP1 ) ) THEN
                  ITMP1 = J
               END IF
  200       CONTINUE
            TMP1 = WORK( M-I+1 )
            WORK( M-I+1 ) = WORK( ITMP1 )
            WORK( ITMP1 ) = TMP1
            IWORK( IWORK( M+ITMP1 ) ) = M - I + 1
            IWORK( IWORK( 2*M-I+1 ) ) = ITMP1
            ITMP2 = IWORK( 2*M-I+1 )
            IWORK( 2*M-I+1 ) = IWORK( M+ITMP1 )
            IWORK( M+ITMP1 ) = ITMP2
*           IWORK( ITMP1 ) = 1
  210    CONTINUE
         J = 0
         DO 220 I = 1, M
            IF( IWORK( I ).GT.LEXTRA .AND. IWORK( I ).LE.M-REXTRA ) THEN
               J = J + 1
               W( J ) = WORK( IWORK( I ) )
               IBLOCK( J ) = IBLOCK( I )
            END IF
  220    CONTINUE
         M = M - LEXTRA - REXTRA
      END IF
      IF( M.NE.ILAST-IFRST+1 ) THEN
         INFO = 2
      END IF
*
  230 CONTINUE
      CALL BLACS_FREEBUFF( ONEDCONTEXT, 1 )
      CALL BLACS_GRIDEXIT( ONEDCONTEXT )
      RETURN
*
*     End of PDSTEBZ
*
      END
*
      SUBROUTINE PDLAEBZ( IJOB, N, MMAX, MINP, ABSTOL, RELTOL, PIVMIN,
     $                    D, NVAL, INTVL, INTVLCT, MOUT, LSAVE, IEFLAG,
     $                    INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*
*     .. Scalar Arguments ..
      INTEGER            IEFLAG, IJOB, INFO, MINP, MMAX, MOUT, N
      DOUBLE PRECISION   ABSTOL, LSAVE, PIVMIN, RELTOL
*     ..
*     .. Array Arguments ..
      INTEGER            INTVLCT( * ), NVAL( * )
      DOUBLE PRECISION   D( * ), INTVL( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAEBZ contains the iteration loop which computes the eigenvalues
*  contained in the input intervals [ INTVL(2*j-1), INTVL(2*j) ] where
*  j = 1,...,MINP. It uses and computes the function N(w), which is
*  the count of eigenvalues of a symmetric tridiagonal matrix less than
*  or equal to its argument w.
*
*  This is a ScaLAPACK internal subroutine and arguments are not
*  checked for unreasonable values.
*
*  Arguments
*  =========
*
*  IJOB    (input) INTEGER
*          Specifies the computation done by PDLAEBZ
*          = 0 : Find an interval with desired values of N(w) at the
*                endpoints of the interval.
*          = 1 : Find a floating point number contained in the initial
*                interval with a desired value of N(w).
*          = 2 : Perform bisection iteration to find eigenvalues of T.
*
*  N       (input) INTEGER
*          The order of the tridiagonal matrix T. N >= 1.
*
*  MMAX    (input) INTEGER
*          The maximum number of intervals that may be generated. If
*          more than MMAX intervals are generated, then PDLAEBZ will
*          quit with INFO = MMAX+1.
*
*  MINP    (input) INTEGER
*          The initial number of intervals. MINP <= MMAX.
*
*  ABSTOL  (input) DOUBLE PRECISION
*          The minimum (absolute) width of an interval. When an interval
*          is narrower than ABSTOL, or than RELTOL times the larger (in
*          magnitude) endpoint, then it is considered to be sufficiently
*          small, i.e., converged.
*          This must be at least zero.
*
*  RELTOL  (input) DOUBLE PRECISION
*          The minimum relative width of an interval. When an interval
*          is narrower than ABSTOL, or than RELTOL times the larger (in
*          magnitude) endpoint, then it is considered to be sufficiently
*          small, i.e., converged.
*          Note : This should be at least radix*machine epsilon.
*
*  PIVMIN  (input) DOUBLE PRECISION
*          The minimum absolute of a "pivot" in the "paranoid"
*          implementation of the Sturm sequence loop. This must be at
*          least max_j |e(j)^2| *safe_min, and at least safe_min, where
*          safe_min is at least the smallest number that can divide 1.0
*          without overflow.
*          See PDLAPDCT for the "paranoid" implementation of the Sturm
*          sequence loop.
*
*  D       (input) DOUBLE PRECISION array, dimension (2*N - 1)
*          Contains the diagonals and the squares of the off-diagonal
*          elements of the tridiagonal matrix T. These elements are
*          assumed to be interleaved in memory for better cache
*          performance. The diagonal entries of T are in the entries
*          D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
*          entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
*          matrix must be scaled so that its largest entry 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.
*
*  NVAL    (input/output) INTEGER array, dimension (4)
*          If IJOB = 0, the desired values of N(w) are in NVAL(1) and
*                       NVAL(2).
*          If IJOB = 1, NVAL(2) is the desired value of N(w).
*          If IJOB = 2, not referenced.
*          This array will, in general, be reordered on output.
*
*  INTVL   (input/output) DOUBLE PRECISION array, dimension (2*MMAX)
*          The endpoints of the intervals. INTVL(2*j-1) is the left
*          endpoint of the j-th interval, and INTVL(2*j) is the right
*          endpoint of the j-th interval. The input intervals will,
*          in general, be modified, split and reordered by the
*          calculation.
*          On input, INTVL contains the MINP input intervals.
*          On output, INTVL contains the converged intervals.
*
*  INTVLCT (input/output) INTEGER array, dimension (2*MMAX)
*          The counts at the endpoints of the intervals. INTVLCT(2*j-1)
*          is the count at the left endpoint of the j-th interval, i.e.,
*          the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
*          count at the right endpoint of the j-th interval.
*          On input, INTVLCT contains the counts at the endpoints of
*          the MINP input intervals.
*          On output, INTVLCT contains the counts at the endpoints of
*          the converged intervals.
*
*  MOUT    (output) INTEGER
*          The number of intervals output.
*
*  LSAVE   (output) DOUBLE PRECISION
*          If IJOB = 0 or 2, not referenced.
*          If IJOB = 1, this is the largest floating point number
*          encountered which has count N(w) = NVAL(1).
*
*  IEFLAG  (input) INTEGER
*          A flag which indicates whether N(w) should be speeded up by
*          exploiting IEEE Arithmetic.
*
*  INFO    (output) INTEGER
*          = 0        : All intervals converged.
*          = 1 - MMAX : The last INFO intervals did not converge.
*          = MMAX + 1 : More than MMAX intervals were generated.
*
*  =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, INT, LOG, MAX, MIN
*     ..
*     .. External Subroutines ..
      EXTERNAL           PDLAECV, PDLAIECTB, PDLAIECTL, PDLAPDCT
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, TWO, HALF
      PARAMETER          ( ZERO = 0.0D+0, TWO = 2.0D+0,
     $                   HALF = 1.0D+0 / TWO )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ITMAX, J, K, KF, KL, KLNEW, L, LCNT, LREQ,
     $                   NALPHA, NBETA, NMID, RCNT, RREQ
      DOUBLE PRECISION   ALPHA, BETA, MID
*     ..
*     .. Executable Statements ..
*
      KF = 1
      KL = MINP + 1
      INFO = 0
      IF( INTVL( 2 )-INTVL( 1 ).LE.ZERO ) THEN
         INFO = MINP
         MOUT = KF
         RETURN
      END IF
      IF( IJOB.EQ.0 ) THEN
*
*        Check if some input intervals have "converged"
*
         CALL PDLAECV( 0, KF, KL, INTVL, INTVLCT, NVAL,
     $                 MAX( ABSTOL, PIVMIN ), RELTOL )
         IF( KF.GE.KL )
     $      GO TO 60
*
*        Compute upper bound on number of iterations needed
*
         ITMAX = INT( ( LOG( INTVL( 2 )-INTVL( 1 )+PIVMIN )-
     $           LOG( PIVMIN ) ) / LOG( TWO ) ) + 2
*
*        Iteration Loop
*
         DO 20 I = 1, ITMAX
            KLNEW = KL
            DO 10 J = KF, KL - 1
               K = 2*J
*
*           Bisect the interval and find the count at that point
*
               MID = HALF*( INTVL( K-1 )+INTVL( K ) )
               IF( IEFLAG.EQ.0 ) THEN
                  CALL PDLAPDCT( MID, N, D, PIVMIN, NMID )
               ELSE IF( IEFLAG.EQ.1 ) THEN
                  CALL PDLAIECTB( MID, N, D, NMID )
               ELSE
                  CALL PDLAIECTL( MID, N, D, NMID )
               END IF
               LREQ = NVAL( K-1 )
               RREQ = NVAL( K )
               IF( KL.EQ.1 )
     $            NMID = MIN( INTVLCT( K ),
     $                   MAX( INTVLCT( K-1 ), NMID ) )
               IF( NMID.LE.NVAL( K-1 ) ) THEN
                  INTVL( K-1 ) = MID
                  INTVLCT( K-1 ) = NMID
               END IF
               IF( NMID.GE.NVAL( K ) ) THEN
                  INTVL( K ) = MID
                  INTVLCT( K ) = NMID
               END IF
               IF( NMID.GT.LREQ .AND. NMID.LT.RREQ ) THEN
                  L = 2*KLNEW
                  INTVL( L-1 ) = MID
                  INTVL( L ) = INTVL( K )
                  INTVLCT( L-1 ) = NVAL( K )
                  INTVLCT( L ) = INTVLCT( K )
                  INTVL( K ) = MID
                  INTVLCT( K ) = NVAL( K-1 )
                  NVAL( L-1 ) = NVAL( K )
                  NVAL( L ) = NVAL( L-1 )
                  NVAL( K ) = NVAL( K-1 )
                  KLNEW = KLNEW + 1
               END IF
   10       CONTINUE
            KL = KLNEW
            CALL PDLAECV( 0, KF, KL, INTVL, INTVLCT, NVAL,
     $                    MAX( ABSTOL, PIVMIN ), RELTOL )
            IF( KF.GE.KL )
     $         GO TO 60
   20    CONTINUE
      ELSE IF( IJOB.EQ.1 ) THEN
         ALPHA = INTVL( 1 )
         BETA = INTVL( 2 )
         NALPHA = INTVLCT( 1 )
         NBETA = INTVLCT( 2 )
         LSAVE = ALPHA
         LREQ = NVAL( 1 )
         RREQ = NVAL( 2 )
   30    CONTINUE
         IF( NBETA.NE.RREQ .AND. BETA-ALPHA.GT.
     $       MAX( ABSTOL, RELTOL*MAX( ABS( ALPHA ), ABS( BETA ) ) ) )
     $        THEN
*
*           Bisect the interval and find the count at that point
*
            MID = HALF*( ALPHA+BETA )
            IF( IEFLAG.EQ.0 ) THEN
               CALL PDLAPDCT( MID, N, D, PIVMIN, NMID )
            ELSE IF( IEFLAG.EQ.1 ) THEN
               CALL PDLAIECTB( MID, N, D, NMID )
            ELSE
               CALL PDLAIECTL( MID, N, D, NMID )
            END IF
            NMID = MIN( NBETA, MAX( NALPHA, NMID ) )
            IF( NMID.GE.RREQ ) THEN
               BETA = MID
               NBETA = NMID
            ELSE
               ALPHA = MID
               NALPHA = NMID
               IF( NMID.EQ.LREQ )
     $            LSAVE = ALPHA
            END IF
            GO TO 30
         END IF
         KL = KF
         INTVL( 1 ) = ALPHA
         INTVL( 2 ) = BETA
         INTVLCT( 1 ) = NALPHA
         INTVLCT( 2 ) = NBETA
      ELSE IF( IJOB.EQ.2 ) THEN
*
*        Check if some input intervals have "converged"
*
         CALL PDLAECV( 1, KF, KL, INTVL, INTVLCT, NVAL,
     $                 MAX( ABSTOL, PIVMIN ), RELTOL )
         IF( KF.GE.KL )
     $      GO TO 60
*
*        Compute upper bound on number of iterations needed
*
         ITMAX = INT( ( LOG( INTVL( 2 )-INTVL( 1 )+PIVMIN )-
     $           LOG( PIVMIN ) ) / LOG( TWO ) ) + 2
*
*        Iteration Loop
*
         DO 50 I = 1, ITMAX
            KLNEW = KL
            DO 40 J = KF, KL - 1
               K = 2*J
               MID = HALF*( INTVL( K-1 )+INTVL( K ) )
               IF( IEFLAG.EQ.0 ) THEN
                  CALL PDLAPDCT( MID, N, D, PIVMIN, NMID )
               ELSE IF( IEFLAG.EQ.1 ) THEN
                  CALL PDLAIECTB( MID, N, D, NMID )
               ELSE
                  CALL PDLAIECTL( MID, N, D, NMID )
               END IF
               LCNT = INTVLCT( K-1 )
               RCNT = INTVLCT( K )
               NMID = MIN( RCNT, MAX( LCNT, NMID ) )
*
*              Form New Interval(s)
*
               IF( NMID.EQ.LCNT ) THEN
                  INTVL( K-1 ) = MID
               ELSE IF( NMID.EQ.RCNT ) THEN
                  INTVL( K ) = MID
               ELSE IF( KLNEW.LT.MMAX+1 ) THEN
                  L = 2*KLNEW
                  INTVL( L-1 ) = MID
                  INTVL( L ) = INTVL( K )
                  INTVLCT( L-1 ) = NMID
                  INTVLCT( L ) = INTVLCT( K )
                  INTVL( K ) = MID
                  INTVLCT( K ) = NMID
                  KLNEW = KLNEW + 1
               ELSE
                  INFO = MMAX + 1
                  RETURN
               END IF
   40       CONTINUE
            KL = KLNEW
            CALL PDLAECV( 1, KF, KL, INTVL, INTVLCT, NVAL,
     $                    MAX( ABSTOL, PIVMIN ), RELTOL )
            IF( KF.GE.KL )
     $         GO TO 60
   50    CONTINUE
      END IF
   60 CONTINUE
      INFO = MAX( KL-KF, 0 )
      MOUT = KL - 1
      RETURN
*
*     End of PDLAEBZ
*
      END
*
*
      SUBROUTINE PDLAECV( IJOB, KF, KL, INTVL, INTVLCT, NVAL, ABSTOL,
     $                    RELTOL )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*
*     .. Scalar Arguments ..
      INTEGER            IJOB, KF, KL
      DOUBLE PRECISION   ABSTOL, RELTOL
*     ..
*     .. Array Arguments ..
      INTEGER            INTVLCT( * ), NVAL( * )
      DOUBLE PRECISION   INTVL( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAECV checks if the input intervals [ INTVL(2*i-1), INTVL(2*i) ],
*  i = KF, ... , KL-1, have "converged".
*  PDLAECV modifies KF to be the index of the last converged interval,
*  i.e., on output, all intervals [ INTVL(2*i-1), INTVL(2*i) ], i < KF,
*  have converged. Note that the input intervals may be reordered by
*  PDLAECV.
*
*  This is a SCALAPACK internal procedure and arguments are not checked
*  for unreasonable values.
*
*  Arguments
*  =========
*
*  IJOB    (input) INTEGER
*          Specifies the criterion for "convergence" of an interval.
*          = 0 : When an interval is narrower than ABSTOL, or than
*                RELTOL times the larger (in magnitude) endpoint, then
*                it is considered to have "converged".
*          = 1 : When an interval is narrower than ABSTOL, or than
*                RELTOL times the larger (in magnitude) endpoint, or if
*                the counts at the endpoints are identical to the counts
*                specified by NVAL ( see NVAL ) then the interval is
*                considered to have "converged".
*
*  KF      (input/output) INTEGER
*          On input, the index of the first input interval is 2*KF-1.
*          On output, the index of the last converged interval
*          is 2*KF-3.
*
*  KL      (input) INTEGER
*          The index of the last input interval is 2*KL-3.
*
*  INTVL   (input/output) DOUBLE PRECISION array, dimension (2*(KL-KF))
*          The endpoints of the intervals. INTVL(2*j-1) is the left
*          oendpoint f the j-th interval, and INTVL(2*j) is the right
*          endpoint of the j-th interval. The input intervals will,
*          in general, be reordered on output.
*          On input, INTVL contains the KL-KF input intervals.
*          On output, INTVL contains the converged intervals, 1 thru'
*          KF-1, and the unconverged intervals, KF thru' KL-1.
*
*  INTVLCT (input/output) INTEGER array, dimension (2*(KL-KF))
*          The counts at the endpoints of the intervals. INTVLCT(2*j-1)
*          is the count at the left endpoint of the j-th interval, i.e.,
*          the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
*          count at the right endpoint of the j-th interval. This array
*          will, in general, be reordered on output.
*          See the comments in PDLAEBZ for more on the function N(w).
*
*  NVAL    (input/output) INTEGER array, dimension (2*(KL-KF))
*          The desired counts, N(w), at the endpoints of the
*          corresponding intervals.  This array will, in general,
*          be reordered on output.
*
*  ABSTOL  (input) DOUBLE PRECISION
*          The minimum (absolute) width of an interval. When an interval
*          is narrower than ABSTOL, or than RELTOL times the larger (in
*          magnitude) endpoint, then it is considered to be sufficiently
*          small, i.e., converged.
*          Note : This must be at least zero.
*
*  RELTOL  (input) DOUBLE PRECISION
*          The minimum relative width of an interval. When an interval
*          is narrower than ABSTOL, or than RELTOL times the larger (in
*          magnitude) endpoint, then it is considered to be sufficiently
*          small, i.e., converged.
*          Note : This should be at least radix*machine epsilon.
*
*  =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX
*     ..
*     .. Local Scalars ..
      LOGICAL            CONDN
      INTEGER            I, ITMP1, ITMP2, J, K, KFNEW
      DOUBLE PRECISION   TMP1, TMP2, TMP3, TMP4
*     ..
*     .. Executable Statements ..
*
      KFNEW = KF
      DO 10 I = KF, KL - 1
         K = 2*I
         TMP3 = INTVL( K-1 )
         TMP4 = INTVL( K )
         TMP1 = ABS( TMP4-TMP3 )
         TMP2 = MAX( ABS( TMP3 ), ABS( TMP4 ) )
         CONDN = TMP1.LT.MAX( ABSTOL, RELTOL*TMP2 )
         IF( IJOB.EQ.0 )
     $      CONDN = CONDN .OR. ( ( INTVLCT( K-1 ).EQ.NVAL( K-1 ) ) .AND.
     $              INTVLCT( K ).EQ.NVAL( K ) )
         IF( CONDN ) THEN
            IF( I.GT.KFNEW ) THEN
*
*              Reorder Intervals
*
               J = 2*KFNEW
               TMP1 = INTVL( K-1 )
               TMP2 = INTVL( K )
               ITMP1 = INTVLCT( K-1 )
               ITMP2 = INTVLCT( K )
               INTVL( K-1 ) = INTVL( J-1 )
               INTVL( K ) = INTVL( J )
               INTVLCT( K-1 ) = INTVLCT( J-1 )
               INTVLCT( K ) = INTVLCT( J )
               INTVL( J-1 ) = TMP1
               INTVL( J ) = TMP2
               INTVLCT( J-1 ) = ITMP1
               INTVLCT( J ) = ITMP2
               IF( IJOB.EQ.0 ) THEN
                  ITMP1 = NVAL( K-1 )
                  NVAL( K-1 ) = NVAL( J-1 )
                  NVAL( J-1 ) = ITMP1
                  ITMP1 = NVAL( K )
                  NVAL( K ) = NVAL( J )
                  NVAL( J ) = ITMP1
               END IF
            END IF
            KFNEW = KFNEW + 1
         END IF
   10 CONTINUE
      KF = KFNEW
      RETURN
*
*     End of PDLAECV
*
      END
*
      SUBROUTINE PDLAPDCT( SIGMA, N, D, PIVMIN, COUNT )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*
*     .. Scalar Arguments ..
      INTEGER            COUNT, N
      DOUBLE PRECISION   PIVMIN, SIGMA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   D( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAPDCT counts the number of negative eigenvalues of (T - SIGMA I).
*  This implementation of the Sturm Sequence loop has conditionals in
*  the innermost loop to avoid overflow and determine the sign of a
*  floating point number. PDLAPDCT will be referred to as the "paranoid"
*  implementation of the Sturm Sequence loop.
*
*  This is a SCALAPACK internal procedure and arguments are not checked
*  for unreasonable values.
*
*  Arguments
*  =========
*
*  SIGMA   (input) DOUBLE PRECISION
*          The shift. PDLAPDCT finds the number of eigenvalues of T less
*          than or equal to SIGMA.
*
*  N       (input) INTEGER
*          The order of the tridiagonal matrix T. N >= 1.
*
*  D       (input) DOUBLE PRECISION array, dimension (2*N - 1)
*          Contains the diagonals and the squares of the off-diagonal
*          elements of the tridiagonal matrix T. These elements are
*          assumed to be interleaved in memory for better cache
*          performance. The diagonal entries of T are in the entries
*          D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
*          entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
*          matrix must be scaled so that its largest entry 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.
*
*  PIVMIN  (input) DOUBLE PRECISION
*          The minimum absolute of a "pivot" in this "paranoid"
*          implementation of the Sturm sequence loop. This must be at
*          least max_j |e(j)^2| *safe_min, and at least safe_min, where
*          safe_min is at least the smallest number that can divide 1.0
*          without overflow.
*
*  COUNT   (output) INTEGER
*          The count of the number of eigenvalues of T less than or
*          equal to SIGMA.
*
*  =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          ABS
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I
      DOUBLE PRECISION   TMP
*     ..
*     .. Executable Statements ..
*
      TMP = D( 1 ) - SIGMA
      IF( ABS( TMP ).LE.PIVMIN )
     $   TMP = -PIVMIN
      COUNT = 0
      IF( TMP.LE.ZERO )
     $   COUNT = 1
      DO 10 I = 3, 2*N - 1, 2
         TMP = D( I ) - D( I-1 ) / TMP - SIGMA
         IF( ABS( TMP ).LE.PIVMIN )
     $      TMP = -PIVMIN
         IF( TMP.LE.ZERO )
     $      COUNT = COUNT + 1
   10 CONTINUE
*
      RETURN
*
*     End of PDLAPDCT
*
      END