Move LAPACK trunk into position.

1 files changed, 916 insertions, 0 deletions

diff --git a/SRC/clarrv.f b/SRC/clarrv.f
new file mode 100644
index 00000000..95cce4e5
--- /dev/null
+++ b/SRC/clarrv.f
@@ -0,0 +1,916 @@
+      SUBROUTINE CLARRV( N, VL, VU, D, L, PIVMIN,
+     $                   ISPLIT, M, DOL, DOU, MINRGP,
+     $                   RTOL1, RTOL2, W, WERR, WGAP,
+     $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
+     $                   WORK, IWORK, INFO )
+*
+*  -- LAPACK auxiliary routine (version 3.1.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            DOL, DOU, INFO, LDZ, M, N
+      REAL               MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
+     $                   ISUPPZ( * ), IWORK( * )
+      REAL               D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
+     $                   WGAP( * ), WORK( * )
+      COMPLEX           Z( LDZ, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CLARRV computes the eigenvectors of the tridiagonal matrix
+*  T = L D L^T given L, D and APPROXIMATIONS to the eigenvalues of L D L^T.
+*  The input eigenvalues should have been computed by SLARRE.
+*
+*  Arguments
+*  =========
+*
+*  N       (input) INTEGER
+*          The order of the matrix.  N >= 0.
+*
+*  VL      (input) REAL            
+*  VU      (input) REAL            
+*          Lower and upper bounds of the interval that contains the desired
+*          eigenvalues. VL < VU. Needed to compute gaps on the left or right
+*          end of the extremal eigenvalues in the desired RANGE.
+*
+*  D       (input/output) REAL             array, dimension (N)
+*          On entry, the N diagonal elements of the diagonal matrix D.
+*          On exit, D may be overwritten.
+*
+*  L       (input/output) REAL             array, dimension (N)
+*          On entry, the (N-1) subdiagonal elements of the unit
+*          bidiagonal matrix L are in elements 1 to N-1 of L
+*          (if the matrix is not splitted.) At the end of each block
+*          is stored the corresponding shift as given by SLARRE.
+*          On exit, L is overwritten.
+*
+*  PIVMIN  (in) DOUBLE PRECISION
+*          The minimum pivot allowed in the Sturm sequence.
+*
+*  ISPLIT  (input) 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.
+*
+*  M       (input) INTEGER
+*          The total number of input eigenvalues.  0 <= M <= N.
+*
+*  DOL     (input) INTEGER
+*  DOU     (input) INTEGER
+*          If the user wants to compute only selected eigenvectors from all
+*          the eigenvalues supplied, he can specify an index range DOL:DOU.
+*          Or else the setting DOL=1, DOU=M should be applied.
+*          Note that DOL and DOU refer to the order in which the eigenvalues
+*          are stored in W.
+*          If the user wants to compute only selected eigenpairs, then
+*          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
+*          computed eigenvectors. All other columns of Z are set to zero.
+*
+*  MINRGP  (input) REAL            
+*
+*  RTOL1   (input) REAL            
+*  RTOL2   (input) REAL            
+*           Parameters for bisection.
+*           An interval [LEFT,RIGHT] has converged if
+*           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
+*
+*  W       (input/output) REAL             array, dimension (N)
+*          The first M elements of W contain the APPROXIMATE eigenvalues for
+*          which eigenvectors are to be computed.  The eigenvalues
+*          should be grouped by split-off block and ordered from
+*          smallest to largest within the block ( The output array
+*          W from SLARRE is expected here ). Furthermore, they are with
+*          respect to the shift of the corresponding root representation
+*          for their block. On exit, W holds the eigenvalues of the
+*          UNshifted matrix.
+*
+*  WERR    (input/output) REAL             array, dimension (N)
+*          The first M elements contain the semiwidth of the uncertainty
+*          interval of the corresponding eigenvalue in W
+*
+*  WGAP    (input/output) REAL             array, dimension (N)
+*          The separation from the right neighbor eigenvalue in W.
+*
+*  IBLOCK  (input) 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  (input) 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 the second block.
+*
+*  GERS    (input) REAL             array, dimension (2*N)
+*          The N Gerschgorin intervals (the i-th Gerschgorin interval
+*          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
+*          be computed from the original UNshifted matrix.
+*
+*  Z       (output) COMPLEX          array, dimension (LDZ, max(1,M) )
+*          If INFO = 0, the first M columns of Z contain the
+*          orthonormal eigenvectors of the matrix T
+*          corresponding to the input eigenvalues, with the i-th
+*          column of Z holding the eigenvector associated with W(i).
+*          Note: the user must ensure that at least max(1,M) columns are
+*          supplied in the array Z.
+*
+*  LDZ     (input) INTEGER
+*          The leading dimension of the array Z.  LDZ >= 1, and if
+*          JOBZ = 'V', LDZ >= max(1,N).
+*
+*  ISUPPZ  (output) INTEGER array, dimension ( 2*max(1,M) )
+*          The support of the eigenvectors in Z, i.e., the indices
+*          indicating the nonzero elements in Z. The I-th eigenvector
+*          is nonzero only in elements ISUPPZ( 2*I-1 ) through
+*          ISUPPZ( 2*I ).
+*
+*  WORK    (workspace) REAL             array, dimension (12*N)
+*
+*  IWORK   (workspace) INTEGER array, dimension (7*N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*
+*          > 0:  A problem occured in CLARRV.
+*          < 0:  One of the called subroutines signaled an internal problem.
+*                Needs inspection of the corresponding parameter IINFO
+*                for further information.
+*
+*          =-1:  Problem in SLARRB when refining a child's eigenvalues.
+*          =-2:  Problem in SLARRF when computing the RRR of a child.
+*                When a child is inside a tight cluster, it can be difficult
+*                to find an RRR. A partial remedy from the user's point of
+*                view is to make the parameter MINRGP smaller and recompile.
+*                However, as the orthogonality of the computed vectors is
+*                proportional to 1/MINRGP, the user should be aware that
+*                he might be trading in precision when he decreases MINRGP.
+*          =-3:  Problem in SLARRB when refining a single eigenvalue
+*                after the Rayleigh correction was rejected.
+*          = 5:  The Rayleigh Quotient Iteration failed to converge to
+*                full accuracy in MAXITR steps.
+*
+*  Further Details
+*  ===============
+*
+*  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
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      INTEGER            MAXITR
+      PARAMETER          ( MAXITR = 10 )
+      COMPLEX            CZERO
+      PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ) )
+      REAL               ZERO, ONE, TWO, THREE, FOUR, HALF
+      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0,
+     $                     TWO = 2.0E0, THREE = 3.0E0,
+     $                     FOUR = 4.0E0, HALF = 0.5E0)
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
+      INTEGER            DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
+     $                   IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
+     $                   INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
+     $                   ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
+     $                   NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
+     $                   NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
+     $                   OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
+     $                   WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
+     $                   ZUSEDW
+      INTEGER            INDIN1, INDIN2
+      REAL               BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
+     $                   LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
+     $                   RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
+     $                   SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
+*     ..
+*     .. External Functions ..
+      REAL               SLAMCH
+      EXTERNAL           SLAMCH
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLAR1V, CLASET, CSSCAL, SCOPY, SLARRB,
+     $                   SLARRF
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC ABS, REAL, MAX, MIN
+      INTRINSIC CMPLX
+*     ..
+*     .. Executable Statements ..
+*     ..
+
+*     The first N entries of WORK are reserved for the eigenvalues
+      INDLD = N+1
+      INDLLD= 2*N+1
+      INDIN1 = 3*N + 1
+      INDIN2 = 4*N + 1
+      INDWRK = 5*N + 1
+      MINWSIZE = 12 * N
+
+      DO 5 I= 1,MINWSIZE
+         WORK( I ) = ZERO
+ 5    CONTINUE
+
+*     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
+*     factorization used to compute the FP vector
+      IINDR = 0
+*     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
+*     layer and the one above.
+      IINDC1 = N
+      IINDC2 = 2*N
+      IINDWK = 3*N + 1
+
+      MINIWSIZE = 7 * N
+      DO 10 I= 1,MINIWSIZE
+         IWORK( I ) = 0
+ 10   CONTINUE
+
+      ZUSEDL = 1
+      IF(DOL.GT.1) THEN
+*        Set lower bound for use of Z
+         ZUSEDL = DOL-1
+      ENDIF
+      ZUSEDU = M
+      IF(DOU.LT.M) THEN
+*        Set lower bound for use of Z
+         ZUSEDU = DOU+1
+      ENDIF
+*     The width of the part of Z that is used
+      ZUSEDW = ZUSEDU - ZUSEDL + 1
+
+
+      CALL CLASET( 'Full', N, ZUSEDW, CZERO, CZERO,
+     $                    Z(1,ZUSEDL), LDZ )
+
+      EPS = SLAMCH( 'Precision' )
+      RQTOL = TWO * EPS
+*
+*     Set expert flags for standard code.
+      TRYRQC = .TRUE.
+
+      IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
+      ELSE
+*        Only selected eigenpairs are computed. Since the other evalues
+*        are not refined by RQ iteration, bisection has to compute to full
+*        accuracy.
+         RTOL1 = FOUR * EPS
+         RTOL2 = FOUR * EPS
+      ENDIF
+
+*     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
+*     desired eigenvalues. The support of the nonzero eigenvector
+*     entries is contained in the interval IBEGIN:IEND.
+*     Remark that if k eigenpairs are desired, then the eigenvectors
+*     are stored in k contiguous columns of Z.
+
+*     DONE is the number of eigenvectors already computed
+      DONE = 0
+      IBEGIN = 1
+      WBEGIN = 1
+      DO 170 JBLK = 1, IBLOCK( M )
+         IEND = ISPLIT( JBLK )
+         SIGMA = L( IEND )
+*        Find the eigenvectors of the submatrix indexed IBEGIN
+*        through IEND.
+         WEND = WBEGIN - 1
+ 15      CONTINUE
+         IF( WEND.LT.M ) THEN
+            IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
+               WEND = WEND + 1
+               GO TO 15
+            END IF
+         END IF
+         IF( WEND.LT.WBEGIN ) THEN
+            IBEGIN = IEND + 1
+            GO TO 170
+         ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
+            IBEGIN = IEND + 1
+            WBEGIN = WEND + 1
+            GO TO 170
+         END IF
+
+*        Find local spectral diameter of the block
+         GL = GERS( 2*IBEGIN-1 )
+         GU = GERS( 2*IBEGIN )
+         DO 20 I = IBEGIN+1 , IEND
+            GL = MIN( GERS( 2*I-1 ), GL )
+            GU = MAX( GERS( 2*I ), GU )
+ 20      CONTINUE
+         SPDIAM = GU - GL
+
+*        OLDIEN is the last index of the previous block
+         OLDIEN = IBEGIN - 1
+*        Calculate the size of the current block
+         IN = IEND - IBEGIN + 1
+*        The number of eigenvalues in the current block
+         IM = WEND - WBEGIN + 1
+
+*        This is for a 1x1 block
+         IF( IBEGIN.EQ.IEND ) THEN
+            DONE = DONE+1
+            Z( IBEGIN, WBEGIN ) = CMPLX( ONE, ZERO )
+            ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
+            ISUPPZ( 2*WBEGIN ) = IBEGIN
+            W( WBEGIN ) = W( WBEGIN ) + SIGMA
+            WORK( WBEGIN ) = W( WBEGIN )
+            IBEGIN = IEND + 1
+            WBEGIN = WBEGIN + 1
+            GO TO 170
+         END IF
+
+*        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
+*        Note that these can be approximations, in this case, the corresp.
+*        entries of WERR give the size of the uncertainty interval.
+*        The eigenvalue approximations will be refined when necessary as
+*        high relative accuracy is required for the computation of the
+*        corresponding eigenvectors.
+         CALL SCOPY( IM, W( WBEGIN ), 1,
+     &                   WORK( WBEGIN ), 1 )
+
+*        We store in W the eigenvalue approximations w.r.t. the original
+*        matrix T.
+         DO 30 I=1,IM
+            W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
+ 30      CONTINUE
+
+
+*        NDEPTH is the current depth of the representation tree
+         NDEPTH = 0
+*        PARITY is either 1 or 0
+         PARITY = 1
+*        NCLUS is the number of clusters for the next level of the
+*        representation tree, we start with NCLUS = 1 for the root
+         NCLUS = 1
+         IWORK( IINDC1+1 ) = 1
+         IWORK( IINDC1+2 ) = IM
+
+*        IDONE is the number of eigenvectors already computed in the current
+*        block
+         IDONE = 0
+*        loop while( IDONE.LT.IM )
+*        generate the representation tree for the current block and
+*        compute the eigenvectors
+   40    CONTINUE
+         IF( IDONE.LT.IM ) THEN
+*           This is a crude protection against infinitely deep trees
+            IF( NDEPTH.GT.M ) THEN
+               INFO = -2
+               RETURN
+            ENDIF
+*           breadth first processing of the current level of the representation
+*           tree: OLDNCL = number of clusters on current level
+            OLDNCL = NCLUS
+*           reset NCLUS to count the number of child clusters
+            NCLUS = 0
+*
+            PARITY = 1 - PARITY
+            IF( PARITY.EQ.0 ) THEN
+               OLDCLS = IINDC1
+               NEWCLS = IINDC2
+            ELSE
+               OLDCLS = IINDC2
+               NEWCLS = IINDC1
+            END IF
+*           Process the clusters on the current level
+            DO 150 I = 1, OLDNCL
+               J = OLDCLS + 2*I
+*              OLDFST, OLDLST = first, last index of current cluster.
+*                               cluster indices start with 1 and are relative
+*                               to WBEGIN when accessing W, WGAP, WERR, Z
+               OLDFST = IWORK( J-1 )
+               OLDLST = IWORK( J )
+               IF( NDEPTH.GT.0 ) THEN
+*                 Retrieve relatively robust representation (RRR) of cluster
+*                 that has been computed at the previous level
+*                 The RRR is stored in Z and overwritten once the eigenvectors
+*                 have been computed or when the cluster is refined
+
+                  IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
+*                    Get representation from location of the leftmost evalue
+*                    of the cluster
+                     J = WBEGIN + OLDFST - 1
+                  ELSE
+                     IF(WBEGIN+OLDFST-1.LT.DOL) THEN
+*                       Get representation from the left end of Z array
+                        J = DOL - 1
+                     ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
+*                       Get representation from the right end of Z array
+                        J = DOU
+                     ELSE
+                        J = WBEGIN + OLDFST - 1
+                     ENDIF
+                  ENDIF
+                  DO 45 K = 1, IN - 1
+                     D( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1,
+     $                                 J ) )
+                     L( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1,
+     $                                 J+1 ) )
+   45             CONTINUE
+                  D( IEND ) = REAL( Z( IEND, J ) )
+                  SIGMA = REAL( Z( IEND, J+1 ) )
+
+*                 Set the corresponding entries in Z to zero
+                  CALL CLASET( 'Full', IN, 2, CZERO, CZERO,
+     $                         Z( IBEGIN, J), LDZ )
+               END IF
+
+*              Compute DL and DLL of current RRR
+               DO 50 J = IBEGIN, IEND-1
+                  TMP = D( J )*L( J )
+                  WORK( INDLD-1+J ) = TMP
+                  WORK( INDLLD-1+J ) = TMP*L( J )
+   50          CONTINUE
+
+               IF( NDEPTH.GT.0 ) THEN
+*                 P and Q are index of the first and last eigenvalue to compute
+*                 within the current block
+                  P = INDEXW( WBEGIN-1+OLDFST )
+                  Q = INDEXW( WBEGIN-1+OLDLST )
+*                 Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET
+*                 thru' Q-OFFSET elements of these arrays are to be used.
+C                  OFFSET = P-OLDFST
+                  OFFSET = INDEXW( WBEGIN ) - 1
+*                 perform limited bisection (if necessary) to get approximate
+*                 eigenvalues to the precision needed.
+                  CALL SLARRB( IN, D( IBEGIN ),
+     $                         WORK(INDLLD+IBEGIN-1),
+     $                         P, Q, RTOL1, RTOL2, OFFSET,
+     $                         WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
+     $                         WORK( INDWRK ), IWORK( IINDWK ),
+     $                         PIVMIN, SPDIAM, IN, IINFO )
+                  IF( IINFO.NE.0 ) THEN
+                     INFO = -1
+                     RETURN
+                  ENDIF
+*                 We also recompute the extremal gaps. W holds all eigenvalues
+*                 of the unshifted matrix and must be used for computation
+*                 of WGAP, the entries of WORK might stem from RRRs with
+*                 different shifts. The gaps from WBEGIN-1+OLDFST to
+*                 WBEGIN-1+OLDLST are correctly computed in SLARRB.
+*                 However, we only allow the gaps to become greater since
+*                 this is what should happen when we decrease WERR
+                  IF( OLDFST.GT.1) THEN
+                     WGAP( WBEGIN+OLDFST-2 ) =
+     $             MAX(WGAP(WBEGIN+OLDFST-2),
+     $                 W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
+     $                 - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
+                  ENDIF
+                  IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
+                     WGAP( WBEGIN+OLDLST-1 ) =
+     $               MAX(WGAP(WBEGIN+OLDLST-1),
+     $                   W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
+     $                   - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
+                  ENDIF
+*                 Each time the eigenvalues in WORK get refined, we store
+*                 the newly found approximation with all shifts applied in W
+                  DO 53 J=OLDFST,OLDLST
+                     W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
+ 53               CONTINUE
+               END IF
+
+*              Process the current node.
+               NEWFST = OLDFST
+               DO 140 J = OLDFST, OLDLST
+                  IF( J.EQ.OLDLST ) THEN
+*                    we are at the right end of the cluster, this is also the
+*                    boundary of the child cluster
+                     NEWLST = J
+                  ELSE IF ( WGAP( WBEGIN + J -1).GE.
+     $                    MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
+*                    the right relative gap is big enough, the child cluster
+*                    (NEWFST,..,NEWLST) is well separated from the following
+                     NEWLST = J
+                   ELSE
+*                    inside a child cluster, the relative gap is not
+*                    big enough.
+                     GOTO 140
+                  END IF
+
+*                 Compute size of child cluster found
+                  NEWSIZ = NEWLST - NEWFST + 1
+
+*                 NEWFTT is the place in Z where the new RRR or the computed
+*                 eigenvector is to be stored
+                  IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
+*                    Store representation at location of the leftmost evalue
+*                    of the cluster
+                     NEWFTT = WBEGIN + NEWFST - 1
+                  ELSE
+                     IF(WBEGIN+NEWFST-1.LT.DOL) THEN
+*                       Store representation at the left end of Z array
+                        NEWFTT = DOL - 1
+                     ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
+*                       Store representation at the right end of Z array
+                        NEWFTT = DOU
+                     ELSE
+                        NEWFTT = WBEGIN + NEWFST - 1
+                     ENDIF
+                  ENDIF
+
+                  IF( NEWSIZ.GT.1) THEN
+*
+*                    Current child is not a singleton but a cluster.
+*                    Compute and store new representation of child.
+*
+*
+*                    Compute left and right cluster gap.
+*
+*                    LGAP and RGAP are not computed from WORK because
+*                    the eigenvalue approximations may stem from RRRs
+*                    different shifts. However, W hold all eigenvalues
+*                    of the unshifted matrix. Still, the entries in WGAP
+*                    have to be computed from WORK since the entries
+*                    in W might be of the same order so that gaps are not
+*                    exhibited correctly for very close eigenvalues.
+                     IF( NEWFST.EQ.1 ) THEN
+                        LGAP = MAX( ZERO,
+     $                       W(WBEGIN)-WERR(WBEGIN) - VL )
+                    ELSE
+                        LGAP = WGAP( WBEGIN+NEWFST-2 )
+                     ENDIF
+                     RGAP = WGAP( WBEGIN+NEWLST-1 )
+*
+*                    Compute left- and rightmost eigenvalue of child
+*                    to high precision in order to shift as close
+*                    as possible and obtain as large relative gaps
+*                    as possible
+*
+                     DO 55 K =1,2
+                        IF(K.EQ.1) THEN
+                           P = INDEXW( WBEGIN-1+NEWFST )
+                        ELSE
+                           P = INDEXW( WBEGIN-1+NEWLST )
+                        ENDIF
+                        OFFSET = INDEXW( WBEGIN ) - 1
+                        CALL SLARRB( IN, D(IBEGIN),
+     $                       WORK( INDLLD+IBEGIN-1 ),P,P,
+     $                       RQTOL, RQTOL, OFFSET,
+     $                       WORK(WBEGIN),WGAP(WBEGIN),
+     $                       WERR(WBEGIN),WORK( INDWRK ),
+     $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
+     $                       IN, IINFO )
+ 55                  CONTINUE
+*
+                     IF((WBEGIN+NEWLST-1.LT.DOL).OR.
+     $                  (WBEGIN+NEWFST-1.GT.DOU)) THEN
+*                       if the cluster contains no desired eigenvalues
+*                       skip the computation of that branch of the rep. tree
+*
+*                       We could skip before the refinement of the extremal
+*                       eigenvalues of the child, but then the representation
+*                       tree could be different from the one when nothing is
+*                       skipped. For this reason we skip at this place.
+                        IDONE = IDONE + NEWLST - NEWFST + 1
+                        GOTO 139
+                     ENDIF
+*
+*                    Compute RRR of child cluster.
+*                    Note that the new RRR is stored in Z
+*
+C                    SLARRF needs LWORK = 2*N
+                     CALL SLARRF( IN, D( IBEGIN ), L( IBEGIN ),
+     $                         WORK(INDLD+IBEGIN-1),
+     $                         NEWFST, NEWLST, WORK(WBEGIN),
+     $                         WGAP(WBEGIN), WERR(WBEGIN),
+     $                         SPDIAM, LGAP, RGAP, PIVMIN, TAU,
+     $                         WORK( INDIN1 ), WORK( INDIN2 ),
+     $                         WORK( INDWRK ), IINFO )
+*                    In the complex case, SLARRF cannot write
+*                    the new RRR directly into Z and needs an intermediate
+*                    workspace
+                     DO 56 K = 1, IN-1
+                        Z( IBEGIN+K-1, NEWFTT ) =
+     $                     CMPLX( WORK( INDIN1+K-1 ), ZERO )
+                        Z( IBEGIN+K-1, NEWFTT+1 ) =
+     $                     CMPLX( WORK( INDIN2+K-1 ), ZERO )
+   56                CONTINUE
+                     Z( IEND, NEWFTT ) =
+     $                  CMPLX( WORK( INDIN1+IN-1 ), ZERO )
+                     IF( IINFO.EQ.0 ) THEN
+*                       a new RRR for the cluster was found by SLARRF
+*                       update shift and store it
+                        SSIGMA = SIGMA + TAU
+                        Z( IEND, NEWFTT+1 ) = CMPLX( SSIGMA, ZERO )
+*                       WORK() are the midpoints and WERR() the semi-width
+*                       Note that the entries in W are unchanged.
+                        DO 116 K = NEWFST, NEWLST
+                           FUDGE =
+     $                          THREE*EPS*ABS(WORK(WBEGIN+K-1))
+                           WORK( WBEGIN + K - 1 ) =
+     $                          WORK( WBEGIN + K - 1) - TAU
+                           FUDGE = FUDGE +
+     $                          FOUR*EPS*ABS(WORK(WBEGIN+K-1))
+*                          Fudge errors
+                           WERR( WBEGIN + K - 1 ) =
+     $                          WERR( WBEGIN + K - 1 ) + FUDGE
+*                          Gaps are not fudged. Provided that WERR is small
+*                          when eigenvalues are close, a zero gap indicates
+*                          that a new representation is needed for resolving
+*                          the cluster. A fudge could lead to a wrong decision
+*                          of judging eigenvalues 'separated' which in
+*                          reality are not. This could have a negative impact
+*                          on the orthogonality of the computed eigenvectors.
+ 116                    CONTINUE
+
+                        NCLUS = NCLUS + 1
+                        K = NEWCLS + 2*NCLUS
+                        IWORK( K-1 ) = NEWFST
+                        IWORK( K ) = NEWLST
+                     ELSE
+                        INFO = -2
+                        RETURN
+                     ENDIF
+                  ELSE
+*
+*                    Compute eigenvector of singleton
+*
+                     ITER = 0
+*
+                     TOL = FOUR * LOG(REAL(IN)) * EPS
+*
+                     K = NEWFST
+                     WINDEX = WBEGIN + K - 1
+                     WINDMN = MAX(WINDEX - 1,1)
+                     WINDPL = MIN(WINDEX + 1,M)
+                     LAMBDA = WORK( WINDEX )
+                     DONE = DONE + 1
+*                    Check if eigenvector computation is to be skipped
+                     IF((WINDEX.LT.DOL).OR.
+     $                  (WINDEX.GT.DOU)) THEN
+                        ESKIP = .TRUE.
+                        GOTO 125
+                     ELSE
+                        ESKIP = .FALSE.
+                     ENDIF
+                     LEFT = WORK( WINDEX ) - WERR( WINDEX )
+                     RIGHT = WORK( WINDEX ) + WERR( WINDEX )
+                     INDEIG = INDEXW( WINDEX )
+*                    Note that since we compute the eigenpairs for a child,
+*                    all eigenvalue approximations are w.r.t the same shift.
+*                    In this case, the entries in WORK should be used for
+*                    computing the gaps since they exhibit even very small
+*                    differences in the eigenvalues, as opposed to the
+*                    entries in W which might "look" the same.
+
+                     IF( K .EQ. 1) THEN
+*                       In the case RANGE='I' and with not much initial
+*                       accuracy in LAMBDA and VL, the formula
+*                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
+*                       can lead to an overestimation of the left gap and
+*                       thus to inadequately early RQI 'convergence'.
+*                       Prevent this by forcing a small left gap.
+                        LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
+                     ELSE
+                        LGAP = WGAP(WINDMN)
+                     ENDIF
+                     IF( K .EQ. IM) THEN
+*                       In the case RANGE='I' and with not much initial
+*                       accuracy in LAMBDA and VU, the formula
+*                       can lead to an overestimation of the right gap and
+*                       thus to inadequately early RQI 'convergence'.
+*                       Prevent this by forcing a small right gap.
+                        RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
+                     ELSE
+                        RGAP = WGAP(WINDEX)
+                     ENDIF
+                     GAP = MIN( LGAP, RGAP )
+                     IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN
+*                       The eigenvector support can become wrong
+*                       because significant entries could be cut off due to a
+*                       large GAPTOL parameter in LAR1V. Prevent this.
+                        GAPTOL = ZERO
+                     ELSE
+                        GAPTOL = GAP * EPS
+                     ENDIF
+                     ISUPMN = IN
+                     ISUPMX = 1
+*                    Update WGAP so that it holds the minimum gap
+*                    to the left or the right. This is crucial in the
+*                    case where bisection is used to ensure that the
+*                    eigenvalue is refined up to the required precision.
+*                    The correct value is restored afterwards.
+                     SAVGAP = WGAP(WINDEX)
+                     WGAP(WINDEX) = GAP
+*                    We want to use the Rayleigh Quotient Correction
+*                    as often as possible since it converges quadratically
+*                    when we are close enough to the desired eigenvalue.
+*                    However, the Rayleigh Quotient can have the wrong sign
+*                    and lead us away from the desired eigenvalue. In this
+*                    case, the best we can do is to use bisection.
+                     USEDBS = .FALSE.
+                     USEDRQ = .FALSE.
+*                    Bisection is initially turned off unless it is forced
+                     NEEDBS =  .NOT.TRYRQC
+ 120                 CONTINUE
+*                    Check if bisection should be used to refine eigenvalue
+                     IF(NEEDBS) THEN
+*                       Take the bisection as new iterate
+                        USEDBS = .TRUE.
+                        ITMP1 = IWORK( IINDR+WINDEX )
+                        OFFSET = INDEXW( WBEGIN ) - 1
+                        CALL SLARRB( IN, D(IBEGIN),
+     $                       WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
+     $                       ZERO, TWO*EPS, OFFSET,
+     $                       WORK(WBEGIN),WGAP(WBEGIN),
+     $                       WERR(WBEGIN),WORK( INDWRK ),
+     $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
+     $                       ITMP1, IINFO )
+                        IF( IINFO.NE.0 ) THEN
+                           INFO = -3
+                           RETURN
+                        ENDIF
+                        LAMBDA = WORK( WINDEX )
+*                       Reset twist index from inaccurate LAMBDA to
+*                       force computation of true MINGMA
+                        IWORK( IINDR+WINDEX ) = 0
+                     ENDIF
+*                    Given LAMBDA, compute the eigenvector.
+                     CALL CLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ),
+     $                    L( IBEGIN ), WORK(INDLD+IBEGIN-1),
+     $                    WORK(INDLLD+IBEGIN-1),
+     $                    PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
+     $                    .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
+     $                    IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
+     $                    NRMINV, RESID, RQCORR, WORK( INDWRK ) )
+                     IF(ITER .EQ. 0) THEN
+                        BSTRES = RESID
+                        BSTW = LAMBDA
+                     ELSEIF(RESID.LT.BSTRES) THEN
+                        BSTRES = RESID
+                        BSTW = LAMBDA
+                     ENDIF
+                     ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
+                     ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
+                     ITER = ITER + 1
+
+*                    sin alpha <= |resid|/gap
+*                    Note that both the residual and the gap are
+*                    proportional to the matrix, so ||T|| doesn't play
+*                    a role in the quotient
+
+*
+*                    Convergence test for Rayleigh-Quotient iteration
+*                    (omitted when Bisection has been used)
+*
+                     IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
+     $                    RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
+     $                    THEN
+*                       We need to check that the RQCORR update doesn't
+*                       move the eigenvalue away from the desired one and
+*                       towards a neighbor. -> protection with bisection
+                        IF(INDEIG.LE.NEGCNT) THEN
+*                          The wanted eigenvalue lies to the left
+                           SGNDEF = -ONE
+                        ELSE
+*                          The wanted eigenvalue lies to the right
+                           SGNDEF = ONE
+                        ENDIF
+*                       We only use the RQCORR if it improves the
+*                       the iterate reasonably.
+                        IF( ( RQCORR*SGNDEF.GE.ZERO )
+     $                       .AND.( LAMBDA + RQCORR.LE. RIGHT)
+     $                       .AND.( LAMBDA + RQCORR.GE. LEFT)
+     $                       ) THEN
+                           USEDRQ = .TRUE.
+*                          Store new midpoint of bisection interval in WORK
+                           IF(SGNDEF.EQ.ONE) THEN
+*                             The current LAMBDA is on the left of the true
+*                             eigenvalue
+                              LEFT = LAMBDA
+*                             We prefer to assume that the error estimate
+*                             is correct. We could make the interval not
+*                             as a bracket but to be modified if the RQCORR
+*                             chooses to. In this case, the RIGHT side should
+*                             be modified as follows:
+*                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
+                           ELSE
+*                             The current LAMBDA is on the right of the true
+*                             eigenvalue
+                              RIGHT = LAMBDA
+*                             See comment about assuming the error estimate is
+*                             correct above.
+*                              LEFT = MIN(LEFT, LAMBDA + RQCORR)
+                           ENDIF
+                           WORK( WINDEX ) =
+     $                       HALF * (RIGHT + LEFT)
+*                          Take RQCORR since it has the correct sign and
+*                          improves the iterate reasonably
+                           LAMBDA = LAMBDA + RQCORR
+*                          Update width of error interval
+                           WERR( WINDEX ) =
+     $                             HALF * (RIGHT-LEFT)
+                        ELSE
+                           NEEDBS = .TRUE.
+                        ENDIF
+                        IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
+*                             The eigenvalue is computed to bisection accuracy
+*                             compute eigenvector and stop
+                           USEDBS = .TRUE.
+                           GOTO 120
+                        ELSEIF( ITER.LT.MAXITR ) THEN
+                           GOTO 120
+                        ELSEIF( ITER.EQ.MAXITR ) THEN
+                           NEEDBS = .TRUE.
+                           GOTO 120
+                        ELSE
+                           INFO = 5
+                           RETURN
+                        END IF
+                     ELSE
+                        STP2II = .FALSE.
+        IF(USEDRQ .AND. USEDBS .AND.
+     $                     BSTRES.LE.RESID) THEN
+                           LAMBDA = BSTW
+                           STP2II = .TRUE.
+                        ENDIF
+                        IF (STP2II) THEN
+*                          improve error angle by second step
+                           CALL CLAR1V( IN, 1, IN, LAMBDA,
+     $                          D( IBEGIN ), L( IBEGIN ),
+     $                          WORK(INDLD+IBEGIN-1),
+     $                          WORK(INDLLD+IBEGIN-1),
+     $                          PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
+     $                          .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
+     $                          IWORK( IINDR+WINDEX ),
+     $                          ISUPPZ( 2*WINDEX-1 ),
+     $                          NRMINV, RESID, RQCORR, WORK( INDWRK ) )
+                        ENDIF
+                        WORK( WINDEX ) = LAMBDA
+                     END IF
+*
+*                    Compute FP-vector support w.r.t. whole matrix
+*
+                     ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
+                     ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
+                     ZFROM = ISUPPZ( 2*WINDEX-1 )
+                     ZTO = ISUPPZ( 2*WINDEX )
+                     ISUPMN = ISUPMN + OLDIEN
+                     ISUPMX = ISUPMX + OLDIEN
+*                    Ensure vector is ok if support in the RQI has changed
+                     IF(ISUPMN.LT.ZFROM) THEN
+                        DO 122 II = ISUPMN,ZFROM-1
+                           Z( II, WINDEX ) = ZERO
+ 122                    CONTINUE
+                     ENDIF
+                     IF(ISUPMX.GT.ZTO) THEN
+                        DO 123 II = ZTO+1,ISUPMX
+                           Z( II, WINDEX ) = ZERO
+ 123                    CONTINUE
+                     ENDIF
+                     CALL CSSCAL( ZTO-ZFROM+1, NRMINV,
+     $                       Z( ZFROM, WINDEX ), 1 )
+ 125                 CONTINUE
+*                    Update W
+                     W( WINDEX ) = LAMBDA+SIGMA
+*                    Recompute the gaps on the left and right
+*                    But only allow them to become larger and not
+*                    smaller (which can only happen through "bad"
+*                    cancellation and doesn't reflect the theory
+*                    where the initial gaps are underestimated due
+*                    to WERR being too crude.)
+                     IF(.NOT.ESKIP) THEN
+                        IF( K.GT.1) THEN
+                           WGAP( WINDMN ) = MAX( WGAP(WINDMN),
+     $                          W(WINDEX)-WERR(WINDEX)
+     $                          - W(WINDMN)-WERR(WINDMN) )
+                        ENDIF
+                        IF( WINDEX.LT.WEND ) THEN
+                           WGAP( WINDEX ) = MAX( SAVGAP,
+     $                          W( WINDPL )-WERR( WINDPL )
+     $                          - W( WINDEX )-WERR( WINDEX) )
+                        ENDIF
+                     ENDIF
+                     IDONE = IDONE + 1
+                  ENDIF
+*                 here ends the code for the current child
+*
+ 139              CONTINUE
+*                 Proceed to any remaining child nodes
+                  NEWFST = J + 1
+ 140           CONTINUE
+ 150        CONTINUE
+            NDEPTH = NDEPTH + 1
+            GO TO 40
+         END IF
+         IBEGIN = IEND + 1
+         WBEGIN = WEND + 1
+ 170  CONTINUE
+*
+
+      RETURN
+*
+*     End of CLARRV
+*
+      END