anchor.diff: +*> \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 )
  -      IMPLICIT NONE
  -*
  -*  -- LAPACK auxiliary routine (version 3.2.2) --
  +*
  +*  -- 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..--
  -*     June 2010
  +*     September 2012
   *
   *     .. Scalar Arguments ..
         CHARACTER          RANGE
  
positive.code:       INTEGER            ISEED( 4 )
  *     ..
  *     .. External Functions ..
        LOGICAL            OLSAME
        DOUBLE PRECISION            ODMCH
        EXTERNAL           ODMCH, OLSAME
  
  *     ..
  *     .. External Subroutines ..
        EXTERNAL           ODCPY, ODRNV, ODRRA, ODRRB, ODRRC, ODRRD,
       $                   ODSQ2
  *     ..
  *     .. Intrinsic Functions ..
  
negative.code:                ISEED( I ) = 1
   122        CONTINUE
  
              CALL DLARNV(2, ISEED, 2*IN-1, WORK(1))
              DO 125 I = 1,IN-1
                 D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
                 E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
  