anchor.diff:        DOUBLE PRECISION   D( * ), E( * ), E2( * ), GERS( * ),
        $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
   *     ..
  -*
  -*  Purpose
  -*  =======
  -*
  -*  To find the desired eigenvalues of a given real symmetric
  -*  tridiagonal matrix T, ODRRE sets any "small" off-diagonal
  -*  elements to zero, and for each unreduced block T_i, it finds
  -*  (a) a suitable shift at one end of the block's spectrum,
  -*  (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
  -*  (c) eigenvalues of each L_i D_i L_i^T.
  -*  The representations and eigenvalues found are then used by
  -*  ODSTMR to compute the eigenvectors of T.
  -*  The accuracy varies depending on whether bisection is used to
  -*  find a few eigenvalues or the dqds algorithm (subroutine ODSQ2) to
  -*  conpute all and then discard any unwanted one.
  -*  As an added benefit, ODRRE also outputs the n
  -*  Gerschgorin intervals for the matrices L_i D_i L_i^T.
  -*
  -*  Arguments
  -*  =========
  -*
  -*  RANGE   (input) CHARACTER
  -*          = 'A': ("All")   all eigenvalues will be found.
  -*          = 'V': ("Value") all eigenvalues in the half-open interval
  -*                           (VL, VU] will be found.
  -*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
  -*                           entire matrix) will be found.
  -*
  -*  N       (input) INTEGER
  -*          The order of the matrix. N > 0.
  -*
  -*  VL      (input/output) DOUBLE PRECISION
  -*  VU      (input/output) DOUBLE PRECISION
  -*          If RANGE='V', the lower and upper bounds for the eigenvalues.
  -*          Eigenvalues less than or equal to VL, or greater than VU,
  -*          will not be returned.  VL < VU.
  -*          If RANGE='I' or ='A', ODRRE computes bounds on the desired
  -*          part of the spectrum.
  -*
  -*  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.
  -*
  -*  D       (input/output) DOUBLE PRECISION array, dimension (N)
  -*          On entry, the N diagonal elements of the tridiagonal
  -*          matrix T.
  -*          On exit, the N diagonal elements of the diagonal
  -*          matrices D_i.
  -*
  -*  E       (input/output) DOUBLE PRECISION array, dimension (N)
  -*          On entry, the first (N-1) entries contain the subdiagonal
  -*          elements of the tridiagonal matrix T; E(N) need not be set.
  -*          On exit, E contains the subdiagonal elements of the unit
  -*          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
  -*          1 <= I <= NSPLIT, contain the base points sigma_i on output.
  -*
  -*  E2      (input/output) DOUBLE PRECISION array, dimension (N)
  -*          On entry, the first (N-1) entries contain the SQUARES of the
  -*          subdiagonal elements of the tridiagonal matrix T;
  -*          E2(N) need not be set.
  -*          On exit, the entries E2( ISPLIT( I ) ),
  -*          1 <= I <= NSPLIT, have been set to zero
  -*
  -*  RTOL1   (input) DOUBLE PRECISION
  -*  RTOL2   (input) DOUBLE PRECISION
  -*           Parameters for bisection.
  -*           An interval [LEFT,RIGHT] has converged if
  -*           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
  -*
  -*  SPLTOL  (input) DOUBLE PRECISION
  -*          The threshold for splitting.
  -*
  -*  NSPLIT  (output) INTEGER
  -*          The number of blocks T splits into. 1 <= NSPLIT <= N.
  -*
  -*  ISPLIT  (output) INTEGER array, dimension (N)
  -*          The splitting points, at which T breaks up into blocks.
  -*          The first block 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.
  -*
  -*  M       (output) INTEGER
  -*          The total number of eigenvalues (of all L_i D_i L_i^T)
  -*          found.
  -*
  -*  W       (output) DOUBLE PRECISION array, dimension (N)
  -*          The first M elements contain the eigenvalues. The
  -*          eigenvalues of each of the blocks, L_i D_i L_i^T, are
  -*          sorted in ascending order ( ODRRE may use the
  -*          remaining N-M elements as workspace).
  -*
  -*  WERR    (output) DOUBLE PRECISION array, dimension (N)
  -*          The error bound on the corresponding eigenvalue in W.
  -*
  -*  WGAP    (output) DOUBLE PRECISION array, dimension (N)
  -*          The separation from the right neighbor eigenvalue in W.
  -*          The gap is only with respect to the eigenvalues of the same block
  -*          as each block has its own representation tree.
  -*          Exception: at the right end of a block we store the left gap
  -*
  -*  IBLOCK  (output) INTEGER array, dimension (N)
  -*          The indices of the blocks (submatrices) associated with the
  -*          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
  -*          W(i) belongs to the first block from the top, =2 if W(i)
  -*          belongs to the second block, etc.
  -*
  -*  INDEXW  (output) INTEGER array, dimension (N)
  -*          The indices of the eigenvalues within each block (submatrix);
  -*          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
  -*          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
  -*
  -*  GERS    (output) DOUBLE PRECISION array, dimension (2*N)
  -*          The N Gerschgorin intervals (the i-th Gerschgorin interval
  -*          is (GERS(2*i-1), GERS(2*i)).
  -*
  -*  PIVMIN  (output) DOUBLE PRECISION
  -*          The minimum pivot in the Sturm sequence for T.
  -*
  -*  WORK    (workspace) DOUBLE PRECISION array, dimension (6*N)
  -*          Workspace.
  -*
  -*  IWORK   (workspace) INTEGER array, dimension (5*N)
  -*          Workspace.
  -*
  -*  INFO    (output) INTEGER
  -*          = 0:  successful exit
  -*          > 0:  A problem occured in ODRRE.
  -*          < 0:  One of the called subroutines signaled an internal problem.
  -*                Needs inspection of the corresponding parameter IINFO
  -*                for further information.
  -*
  -*          =-1:  Problem in ODRRD.
  -*          = 2:  No base representation could be found in MAXTRY iterations.
  -*                Increasing MAXTRY and recompilation might be a remedy.
  -*          =-3:  Problem in ODRRB when computing the refined root
  -*                representation for ODSQ2.
  -*          =-4:  Problem in ODRRB when preforming bisection on the
  -*                desired part of the spectrum.
  -*          =-5:  Problem in ODSQ2.
  -*          =-6:  Problem in ODSQ2.
  -*
  -*  Further Details
  -*  The base representations are required to suffer very little
  -*  element growth and consequently define all their eigenvalues to
  -*  high relative accuracy.
  -*  ===============
  -*
  -*  Based on contributions by
  -*     Beresford Parlett, University of California, Berkeley, USA
  -*     Jim Demmel, University of California, Berkeley, USA
  -*     Inderjit Dhillon, University of Texas, Austin, USA
  -*     Osni Marques, LBNL/NERSC, USA
  -*     Christof Voemel, University of California, Berkeley, USA
   *
   *  =====================================================================
   *
  
positive.code: *        Compute the negcount at the 1/4 and 3/4 points
           IF(MB.GT.1) THEN
              CALL ODRRC( 'T', IN, S1, S2, D(IBEGIN),
       $                   E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
           ENDIF
  
           IF(MB.EQ.1) THEN
  
negative.code: *> \brief \b ODRRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.
  *
  *  =========== DOCUMENTATION ===========
  *
  * Online html documentation available at 
  *            http://www.netlib.org/lapack/explore-html/ 
  *
  *> \htmlonly
  *> Download ODRRE + dependencies 
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f"> 
  *> [TGZ]</a> 
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f"> 
  *> [ZIP]</a> 
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f"> 
  *> [TXT]</a>
  *> \endhtmlonly 
  *
  *  Definition:
  *  ===========
  *
  *       SUBROUTINE ODRRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
  *                           RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
  *                           W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
  *                           WORK, IWORK, INFO )
  * 
  *       .. Scalar Arguments ..
  *       CHARACTER          RANGE
  *       INTEGER            IL, INFO, IU, M, N, NSPLIT
  *       DOUBLE PRECISION  PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
  *       ..
  *       .. Array Arguments ..
  *       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
  *      $                   INDEXW( * )
  *       DOUBLE PRECISION   D( * ), E( * ), E2( * ), GERS( * ),
  *      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
  *       ..
  *  
  *
  *> \par Purpose:
  *  =============
  *>
  *> \verbatim
  *>
  *> To find the desired eigenvalues of a given real symmetric
  *> tridiagonal matrix T, ODRRE sets any "small" off-diagonal
  *> elements to zero, and for each unreduced block T_i, it finds
  *> (a) a suitable shift at one end of the block's spectrum,
  *> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
  *> (c) eigenvalues of each L_i D_i L_i^T.
  *> The representations and eigenvalues found are then used by
  *> DSTEMR to compute the eigenvectors of T.
  *> The accuracy varies depending on whether bisection is used to
  *> find a few eigenvalues or the dqds algorithm (subroutine ODSQ2) to
  *> conpute all and then discard any unwanted one.
  *> As an added benefit, ODRRE also outputs the n
  *> Gerschgorin intervals for the matrices L_i D_i L_i^T.
  *> \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] N
  *> \verbatim
  *>          N is INTEGER
  *>          The order of the matrix. N > 0.
  *> \endverbatim
  *>
  *> \param[in,out] VL
  *> \verbatim
  *>          VL is DOUBLE PRECISION
  *> \endverbatim
  *>
  *> \param[in,out] VU
  *> \verbatim
  *>          VU is DOUBLE PRECISION
  *>          If RANGE='V', the lower and upper bounds for the eigenvalues.
  *>          Eigenvalues less than or equal to VL, or greater than VU,
  *>          will not be returned.  VL < VU.
  *>          If RANGE='I' or ='A', ODRRE computes bounds on the desired
  *>          part of the spectrum.
  *> \endverbatim
  *>
  *> \param[in] IL
  *> \verbatim
  *>          IL is INTEGER
  *> \endverbatim
  *>
  *> \param[in] IU
  *> \verbatim
  *>          IU is INTEGER
  *>          If RANGE='I', the indices (in ascending order) of the
  *>          smallest and largest eigenvalues to be returned.
  *>          1 <= IL <= IU <= N.
  *> \endverbatim
  *>
  *> \param[in,out] D
  *> \verbatim
  *>          D is DOUBLE PRECISION array, dimension (N)
  *>          On entry, the N diagonal elements of the tridiagonal
  *>          matrix T.
  *>          On exit, the N diagonal elements of the diagonal
  *>          matrices D_i.
  *> \endverbatim
  *>
  *> \param[in,out] E
  *> \verbatim
  *>          E is DOUBLE PRECISION array, dimension (N)
  *>          On entry, the first (N-1) entries contain the subdiagonal
  *>          elements of the tridiagonal matrix T; E(N) need not be set.
  *>          On exit, E contains the subdiagonal elements of the unit
  *>          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
  *>          1 <= I <= NSPLIT, contain the base points sigma_i on output.
  *> \endverbatim
  *>
  *> \param[in,out] E2
  *> \verbatim
  *>          E2 is DOUBLE PRECISION array, dimension (N)
  *>          On entry, the first (N-1) entries contain the SQUARES of the
  *>          subdiagonal elements of the tridiagonal matrix T;
  *>          E2(N) need not be set.
  *>          On exit, the entries E2( ISPLIT( I ) ),
  *>          1 <= I <= NSPLIT, have been set to zero
  *> \endverbatim
  *>
  *> \param[in] RTOL1
  *> \verbatim
  *>          RTOL1 is DOUBLE PRECISION
  *> \endverbatim
  *>
  *> \param[in] RTOL2
  *> \verbatim
  *>          RTOL2 is DOUBLE PRECISION
  *>           Parameters for bisection.
  *>           An interval [LEFT,RIGHT] has converged if
  *>           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
  *> \endverbatim
  *>
  *> \param[in] SPLTOL
  *> \verbatim
  *>          SPLTOL is DOUBLE PRECISION
  *>          The threshold for splitting.
  *> \endverbatim
  *>
  *> \param[out] NSPLIT
  *> \verbatim
  *>          NSPLIT is INTEGER
  *>          The number of blocks T splits into. 1 <= NSPLIT <= N.
  *> \endverbatim
  *>
  *> \param[out] ISPLIT
  *> \verbatim
  *>          ISPLIT is INTEGER array, dimension (N)
  *>          The splitting points, at which T breaks up into blocks.
  *>          The first block 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.
  *> \endverbatim
  *>
  *> \param[out] M
  *> \verbatim
  *>          M is INTEGER
  *>          The total number of eigenvalues (of all L_i D_i L_i^T)
  *>          found.
  *> \endverbatim
  *>
  *> \param[out] W
  *> \verbatim
  *>          W is DOUBLE PRECISION array, dimension (N)
  *>          The first M elements contain the eigenvalues. The
  *>          eigenvalues of each of the blocks, L_i D_i L_i^T, are
  *>          sorted in ascending order ( ODRRE may use the
  *>          remaining N-M elements as workspace).
  *> \endverbatim
  *>
  *> \param[out] WERR
  *> \verbatim
  *>          WERR is DOUBLE PRECISION array, dimension (N)
  *>          The error bound on the corresponding eigenvalue in W.
  *> \endverbatim
  *>
  *> \param[out] WGAP
  *> \verbatim
  *>          WGAP is DOUBLE PRECISION array, dimension (N)
  *>          The separation from the right neighbor eigenvalue in W.
  *>          The gap is only with respect to the eigenvalues of the same block
  *>          as each block has its own representation tree.
  *>          Exception: at the right end of a block we store the left gap
  *> \endverbatim
  *>
  *> \param[out] IBLOCK
  *> \verbatim
  *>          IBLOCK is INTEGER array, dimension (N)
  *>          The indices of the blocks (submatrices) associated with the
  *>          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
  *>          W(i) belongs to the first block from the top, =2 if W(i)
  *>          belongs to the second block, etc.
  *> \endverbatim
  *>
  *> \param[out] INDEXW
  *> \verbatim
  *>          INDEXW is INTEGER array, dimension (N)
  *>          The indices of the eigenvalues within each block (submatrix);
  *>          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
  *>          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
  *> \endverbatim
  *>
  *> \param[out] GERS
  *> \verbatim
  *>          GERS is DOUBLE PRECISION array, dimension (2*N)
  *>          The N Gerschgorin intervals (the i-th Gerschgorin interval
  *>          is (GERS(2*i-1), GERS(2*i)).
  *> \endverbatim
  *>
  *> \param[out] PIVMIN
  *> \verbatim
  *>          PIVMIN is DOUBLE PRECISION
  *>          The minimum pivot in the Sturm sequence for T.
  *> \endverbatim
  *>
  *> \param[out] WORK
  *> \verbatim
  *>          WORK is DOUBLE PRECISION array, dimension (6*N)
  *>          Workspace.
  *> \endverbatim
  *>
  *> \param[out] IWORK
  *> \verbatim
  *>          IWORK is INTEGER array, dimension (5*N)
  *>          Workspace.
  *> \endverbatim
  *>
  *> \param[out] INFO
  *> \verbatim
  *>          INFO is INTEGER
  *>          = 0:  successful exit
  *>          > 0:  A problem occured in ODRRE.
  *>          < 0:  One of the called subroutines signaled an internal problem.
  *>                Needs inspection of the corresponding parameter IINFO
  *>                for further information.
  *>
  *>          =-1:  Problem in ODRRD.
  *>          = 2:  No base representation could be found in MAXTRY iterations.
  *>                Increasing MAXTRY and recompilation might be a remedy.
  *>          =-3:  Problem in ODRRB when computing the refined root
  *>                representation for ODSQ2.
  *>          =-4:  Problem in ODRRB when preforming bisection on the
  *>                desired part of the spectrum.
  *>          =-5:  Problem in ODSQ2.
  *>          =-6:  Problem in ODSQ2.
  *> \endverbatim
  *
  *  Authors:
  *  ========
  *
  *> \author Univ. of Tennessee 
  *> \author Univ. of California Berkeley 
  *> \author Univ. of Colorado Denver 
  *> \author NAG Ltd. 
  *
  *> \date September 2012
  *
  *> \ingroup auxOTHERauxiliary
  *
  *> \par Further Details:
  *  =====================
  *>
  *> \verbatim
  *>
  *>  The base representations are required to suffer very little
  *>  element growth and consequently define all their eigenvalues to
  *>  high relative accuracy.
  *> \endverbatim
  *
  *> \par Contributors:
  *  ==================
  *>
  *>     Beresford Parlett, University of California, Berkeley, USA \n
  *>     Jim Demmel, University of California, Berkeley, USA \n
  *>     Inderjit Dhillon, University of Texas, Austin, USA \n
  *>     Osni Marques, LBNL/NERSC, USA \n
  *>     Christof Voemel, University of California, Berkeley, USA \n
  *>
  *  =====================================================================
        SUBROUTINE ODRRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
       $                    RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
       $                    W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
       $                    WORK, IWORK, INFO )
  *
  *  -- LAPACK auxiliary routine (version 3.4.2) --
  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  *     September 2012
  *
  *     .. Scalar Arguments ..
        CHARACTER          RANGE
  