SRC/dlarrv.f

   1 *> \brief \b DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLARRV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
  22 *                          ISPLIT, M, DOL, DOU, MINRGP,
  23 *                          RTOL1, RTOL2, W, WERR, WGAP,
  24 *                          IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
  25 *                          WORK, IWORK, INFO )
  26 *
  27 *       .. Scalar Arguments ..
  28 *       INTEGER            DOL, DOU, INFO, LDZ, M, N
  29 *       DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
  33 *      $                   ISUPPZ( * ), IWORK( * )
  34 *       DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
  35 *      $                   WGAP( * ), WORK( * )
  36 *       DOUBLE PRECISION  Z( LDZ, * )
  37 *       ..
  38 *
  39 *
  40 *> \par Purpose:
  41 *  =============
  42 *>
  43 *> \verbatim
  44 *>
  45 *> DLARRV computes the eigenvectors of the tridiagonal matrix
  46 *> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
  47 *> The input eigenvalues should have been computed by DLARRE.
  48 *> \endverbatim
  49 *
  50 *  Arguments:
  51 *  ==========
  52 *
  53 *> \param[in] N
  54 *> \verbatim
  55 *>          N is INTEGER
  56 *>          The order of the matrix.  N >= 0.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] VL
  60 *> \verbatim
  61 *>          VL is DOUBLE PRECISION
  62 *>          Lower bound of the interval that contains the desired
  63 *>          eigenvalues. VL < VU. Needed to compute gaps on the left or right
  64 *>          end of the extremal eigenvalues in the desired RANGE.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] VU
  68 *> \verbatim
  69 *>          VU is DOUBLE PRECISION
  70 *>          Upper bound of the interval that contains the desired
  71 *>          eigenvalues. VL < VU. Needed to compute gaps on the left or right
  72 *>          end of the extremal eigenvalues in the desired RANGE.
  73 *> \endverbatim
  74 *>
  75 *> \param[in,out] D
  76 *> \verbatim
  77 *>          D is DOUBLE PRECISION array, dimension (N)
  78 *>          On entry, the N diagonal elements of the diagonal matrix D.
  79 *>          On exit, D may be overwritten.
  80 *> \endverbatim
  81 *>
  82 *> \param[in,out] L
  83 *> \verbatim
  84 *>          L is DOUBLE PRECISION array, dimension (N)
  85 *>          On entry, the (N-1) subdiagonal elements of the unit
  86 *>          bidiagonal matrix L are in elements 1 to N-1 of L
  87 *>          (if the matrix is not split.) At the end of each block
  88 *>          is stored the corresponding shift as given by DLARRE.
  89 *>          On exit, L is overwritten.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] PIVMIN
  93 *> \verbatim
  94 *>          PIVMIN is DOUBLE PRECISION
  95 *>          The minimum pivot allowed in the Sturm sequence.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] ISPLIT
  99 *> \verbatim
 100 *>          ISPLIT is INTEGER array, dimension (N)
 101 *>          The splitting points, at which T breaks up into blocks.
 102 *>          The first block consists of rows/columns 1 to
 103 *>          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
 104 *>          through ISPLIT( 2 ), etc.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] M
 108 *> \verbatim
 109 *>          M is INTEGER
 110 *>          The total number of input eigenvalues.  0 <= M <= N.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] DOL
 114 *> \verbatim
 115 *>          DOL is INTEGER
 116 *> \endverbatim
 117 *>
 118 *> \param[in] DOU
 119 *> \verbatim
 120 *>          DOU is INTEGER
 121 *>          If the user wants to compute only selected eigenvectors from all
 122 *>          the eigenvalues supplied, he can specify an index range DOL:DOU.
 123 *>          Or else the setting DOL=1, DOU=M should be applied.
 124 *>          Note that DOL and DOU refer to the order in which the eigenvalues
 125 *>          are stored in W.
 126 *>          If the user wants to compute only selected eigenpairs, then
 127 *>          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
 128 *>          computed eigenvectors. All other columns of Z are set to zero.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] MINRGP
 132 *> \verbatim
 133 *>          MINRGP is DOUBLE PRECISION
 134 *> \endverbatim
 135 *>
 136 *> \param[in] RTOL1
 137 *> \verbatim
 138 *>          RTOL1 is DOUBLE PRECISION
 139 *> \endverbatim
 140 *>
 141 *> \param[in] RTOL2
 142 *> \verbatim
 143 *>          RTOL2 is DOUBLE PRECISION
 144 *>           Parameters for bisection.
 145 *>           An interval [LEFT,RIGHT] has converged if
 146 *>           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
 147 *> \endverbatim
 148 *>
 149 *> \param[in,out] W
 150 *> \verbatim
 151 *>          W is DOUBLE PRECISION array, dimension (N)
 152 *>          The first M elements of W contain the APPROXIMATE eigenvalues for
 153 *>          which eigenvectors are to be computed.  The eigenvalues
 154 *>          should be grouped by split-off block and ordered from
 155 *>          smallest to largest within the block ( The output array
 156 *>          W from DLARRE is expected here ). Furthermore, they are with
 157 *>          respect to the shift of the corresponding root representation
 158 *>          for their block. On exit, W holds the eigenvalues of the
 159 *>          UNshifted matrix.
 160 *> \endverbatim
 161 *>
 162 *> \param[in,out] WERR
 163 *> \verbatim
 164 *>          WERR is DOUBLE PRECISION array, dimension (N)
 165 *>          The first M elements contain the semiwidth of the uncertainty
 166 *>          interval of the corresponding eigenvalue in W
 167 *> \endverbatim
 168 *>
 169 *> \param[in,out] WGAP
 170 *> \verbatim
 171 *>          WGAP is DOUBLE PRECISION array, dimension (N)
 172 *>          The separation from the right neighbor eigenvalue in W.
 173 *> \endverbatim
 174 *>
 175 *> \param[in] IBLOCK
 176 *> \verbatim
 177 *>          IBLOCK is INTEGER array, dimension (N)
 178 *>          The indices of the blocks (submatrices) associated with the
 179 *>          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
 180 *>          W(i) belongs to the first block from the top, =2 if W(i)
 181 *>          belongs to the second block, etc.
 182 *> \endverbatim
 183 *>
 184 *> \param[in] INDEXW
 185 *> \verbatim
 186 *>          INDEXW is INTEGER array, dimension (N)
 187 *>          The indices of the eigenvalues within each block (submatrix);
 188 *>          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
 189 *>          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
 190 *> \endverbatim
 191 *>
 192 *> \param[in] GERS
 193 *> \verbatim
 194 *>          GERS is DOUBLE PRECISION array, dimension (2*N)
 195 *>          The N Gerschgorin intervals (the i-th Gerschgorin interval
 196 *>          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
 197 *>          be computed from the original UNshifted matrix.
 198 *> \endverbatim
 199 *>
 200 *> \param[out] Z
 201 *> \verbatim
 202 *>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
 203 *>          If INFO = 0, the first M columns of Z contain the
 204 *>          orthonormal eigenvectors of the matrix T
 205 *>          corresponding to the input eigenvalues, with the i-th
 206 *>          column of Z holding the eigenvector associated with W(i).
 207 *>          Note: the user must ensure that at least max(1,M) columns are
 208 *>          supplied in the array Z.
 209 *> \endverbatim
 210 *>
 211 *> \param[in] LDZ
 212 *> \verbatim
 213 *>          LDZ is INTEGER
 214 *>          The leading dimension of the array Z.  LDZ >= 1, and if
 215 *>          JOBZ = 'V', LDZ >= max(1,N).
 216 *> \endverbatim
 217 *>
 218 *> \param[out] ISUPPZ
 219 *> \verbatim
 220 *>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
 221 *>          The support of the eigenvectors in Z, i.e., the indices
 222 *>          indicating the nonzero elements in Z. The I-th eigenvector
 223 *>          is nonzero only in elements ISUPPZ( 2*I-1 ) through
 224 *>          ISUPPZ( 2*I ).
 225 *> \endverbatim
 226 *>
 227 *> \param[out] WORK
 228 *> \verbatim
 229 *>          WORK is DOUBLE PRECISION array, dimension (12*N)
 230 *> \endverbatim
 231 *>
 232 *> \param[out] IWORK
 233 *> \verbatim
 234 *>          IWORK is INTEGER array, dimension (7*N)
 235 *> \endverbatim
 236 *>
 237 *> \param[out] INFO
 238 *> \verbatim
 239 *>          INFO is INTEGER
 240 *>          = 0:  successful exit
 241 *>
 242 *>          > 0:  A problem occurred in DLARRV.
 243 *>          < 0:  One of the called subroutines signaled an internal problem.
 244 *>                Needs inspection of the corresponding parameter IINFO
 245 *>                for further information.
 246 *>
 247 *>          =-1:  Problem in DLARRB when refining a child's eigenvalues.
 248 *>          =-2:  Problem in DLARRF when computing the RRR of a child.
 249 *>                When a child is inside a tight cluster, it can be difficult
 250 *>                to find an RRR. A partial remedy from the user's point of
 251 *>                view is to make the parameter MINRGP smaller and recompile.
 252 *>                However, as the orthogonality of the computed vectors is
 253 *>                proportional to 1/MINRGP, the user should be aware that
 254 *>                he might be trading in precision when he decreases MINRGP.
 255 *>          =-3:  Problem in DLARRB when refining a single eigenvalue
 256 *>                after the Rayleigh correction was rejected.
 257 *>          = 5:  The Rayleigh Quotient Iteration failed to converge to
 258 *>                full accuracy in MAXITR steps.
 259 *> \endverbatim
 260 *
 261 *  Authors:
 262 *  ========
 263 *
 264 *> \author Univ. of Tennessee
 265 *> \author Univ. of California Berkeley
 266 *> \author Univ. of Colorado Denver
 267 *> \author NAG Ltd.
 268 *
 269 *> \date June 2016
 270 *
 271 *> \ingroup doubleOTHERauxiliary
 272 *
 273 *> \par Contributors:
 274 *  ==================
 275 *>
 276 *> Beresford Parlett, University of California, Berkeley, USA \n
 277 *> Jim Demmel, University of California, Berkeley, USA \n
 278 *> Inderjit Dhillon, University of Texas, Austin, USA \n
 279 *> Osni Marques, LBNL/NERSC, USA \n
 280 *> Christof Voemel, University of California, Berkeley, USA
 281 *
 282 *  =====================================================================
 283       SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
 284      $                   ISPLIT, M, DOL, DOU, MINRGP,
 285      $                   RTOL1, RTOL2, W, WERR, WGAP,
 286      $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
 287      $                   WORK, IWORK, INFO )
 288 *
 289 *  -- LAPACK auxiliary routine (version 3.6.1) --
 290 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 291 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 292 *     June 2016
 293 *
 294 *     .. Scalar Arguments ..
 295       INTEGER            DOL, DOU, INFO, LDZ, M, N
 296       DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
 297 *     ..
 298 *     .. Array Arguments ..
 299       INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
 300      $                   ISUPPZ( * ), IWORK( * )
 301       DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
 302      $                   WGAP( * ), WORK( * )
 303       DOUBLE PRECISION  Z( LDZ, * )
 304 *     ..
 305 *
 306 *  =====================================================================
 307 *
 308 *     .. Parameters ..
 309       INTEGER            MAXITR
 310       PARAMETER          ( MAXITR = 10 )
 311       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, HALF
 312       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
 313      $                     TWO = 2.0D0, THREE = 3.0D0,
 314      $                     FOUR = 4.0D0, HALF = 0.5D0)
 315 *     ..
 316 *     .. Local Scalars ..
 317       LOGICAL            ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
 318       INTEGER            DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
 319      $                   IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
 320      $                   INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
 321      $                   ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
 322      $                   NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
 323      $                   NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
 324      $                   OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
 325      $                   WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
 326      $                   ZUSEDW
 327       DOUBLE PRECISION   BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
 328      $                   LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
 329      $                   RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
 330      $                   SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
 331 *     ..
 332 *     .. External Functions ..
 333       DOUBLE PRECISION   DLAMCH
 334       EXTERNAL           DLAMCH
 335 *     ..
 336 *     .. External Subroutines ..
 337       EXTERNAL           DCOPY, DLAR1V, DLARRB, DLARRF, DLASET,
 338      $                   DSCAL
 339 *     ..
 340 *     .. Intrinsic Functions ..
 341       INTRINSIC ABS, DBLE, MAX, MIN
 342 *     ..
 343 *     .. Executable Statements ..
 344 *     ..
 345
 346       INFO = 0
 347 *     The first N entries of WORK are reserved for the eigenvalues
 348       INDLD = N+1
 349       INDLLD= 2*N+1
 350       INDWRK= 3*N+1
 351       MINWSIZE = 12 * N
 352
 353       DO 5 I= 1,MINWSIZE
 354          WORK( I ) = ZERO
 355  5    CONTINUE
 356
 357 *     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
 358 *     factorization used to compute the FP vector
 359       IINDR = 0
 360 *     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
 361 *     layer and the one above.
 362       IINDC1 = N
 363       IINDC2 = 2*N
 364       IINDWK = 3*N + 1
 365
 366       MINIWSIZE = 7 * N
 367       DO 10 I= 1,MINIWSIZE
 368          IWORK( I ) = 0
 369  10   CONTINUE
 370
 371       ZUSEDL = 1
 372       IF(DOL.GT.1) THEN
 373 *        Set lower bound for use of Z
 374          ZUSEDL = DOL-1
 375       ENDIF
 376       ZUSEDU = M
 377       IF(DOU.LT.M) THEN
 378 *        Set lower bound for use of Z
 379          ZUSEDU = DOU+1
 380       ENDIF
 381 *     The width of the part of Z that is used
 382       ZUSEDW = ZUSEDU - ZUSEDL + 1
 383
 384
 385       CALL DLASET( 'Full', N, ZUSEDW, ZERO, ZERO,
 386      $                    Z(1,ZUSEDL), LDZ )
 387
 388       EPS = DLAMCH( 'Precision' )
 389       RQTOL = TWO * EPS
 390 *
 391 *     Set expert flags for standard code.
 392       TRYRQC = .TRUE.
 393
 394       IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
 395       ELSE
 396 *        Only selected eigenpairs are computed. Since the other evalues
 397 *        are not refined by RQ iteration, bisection has to compute to full
 398 *        accuracy.
 399          RTOL1 = FOUR * EPS
 400          RTOL2 = FOUR * EPS
 401       ENDIF
 402
 403 *     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
 404 *     desired eigenvalues. The support of the nonzero eigenvector
 405 *     entries is contained in the interval IBEGIN:IEND.
 406 *     Remark that if k eigenpairs are desired, then the eigenvectors
 407 *     are stored in k contiguous columns of Z.
 408
 409 *     DONE is the number of eigenvectors already computed
 410       DONE = 0
 411       IBEGIN = 1
 412       WBEGIN = 1
 413       DO 170 JBLK = 1, IBLOCK( M )
 414          IEND = ISPLIT( JBLK )
 415          SIGMA = L( IEND )
 416 *        Find the eigenvectors of the submatrix indexed IBEGIN
 417 *        through IEND.
 418          WEND = WBEGIN - 1
 419  15      CONTINUE
 420          IF( WEND.LT.M ) THEN
 421             IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
 422                WEND = WEND + 1
 423                GO TO 15
 424             END IF
 425          END IF
 426          IF( WEND.LT.WBEGIN ) THEN
 427             IBEGIN = IEND + 1
 428             GO TO 170
 429          ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
 430             IBEGIN = IEND + 1
 431             WBEGIN = WEND + 1
 432             GO TO 170
 433          END IF
 434
 435 *        Find local spectral diameter of the block
 436          GL = GERS( 2*IBEGIN-1 )
 437          GU = GERS( 2*IBEGIN )
 438          DO 20 I = IBEGIN+1 , IEND
 439             GL = MIN( GERS( 2*I-1 ), GL )
 440             GU = MAX( GERS( 2*I ), GU )
 441  20      CONTINUE
 442          SPDIAM = GU - GL
 443
 444 *        OLDIEN is the last index of the previous block
 445          OLDIEN = IBEGIN - 1
 446 *        Calculate the size of the current block
 447          IN = IEND - IBEGIN + 1
 448 *        The number of eigenvalues in the current block
 449          IM = WEND - WBEGIN + 1
 450
 451 *        This is for a 1x1 block
 452          IF( IBEGIN.EQ.IEND ) THEN
 453             DONE = DONE+1
 454             Z( IBEGIN, WBEGIN ) = ONE
 455             ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
 456             ISUPPZ( 2*WBEGIN ) = IBEGIN
 457             W( WBEGIN ) = W( WBEGIN ) + SIGMA
 458             WORK( WBEGIN ) = W( WBEGIN )
 459             IBEGIN = IEND + 1
 460             WBEGIN = WBEGIN + 1
 461             GO TO 170
 462          END IF
 463
 464 *        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
 465 *        Note that these can be approximations, in this case, the corresp.
 466 *        entries of WERR give the size of the uncertainty interval.
 467 *        The eigenvalue approximations will be refined when necessary as
 468 *        high relative accuracy is required for the computation of the
 469 *        corresponding eigenvectors.
 470          CALL DCOPY( IM, W( WBEGIN ), 1,
 471      $                   WORK( WBEGIN ), 1 )
 472
 473 *        We store in W the eigenvalue approximations w.r.t. the original
 474 *        matrix T.
 475          DO 30 I=1,IM
 476             W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
 477  30      CONTINUE
 478
 479
 480 *        NDEPTH is the current depth of the representation tree
 481          NDEPTH = 0
 482 *        PARITY is either 1 or 0
 483          PARITY = 1
 484 *        NCLUS is the number of clusters for the next level of the
 485 *        representation tree, we start with NCLUS = 1 for the root
 486          NCLUS = 1
 487          IWORK( IINDC1+1 ) = 1
 488          IWORK( IINDC1+2 ) = IM
 489
 490 *        IDONE is the number of eigenvectors already computed in the current
 491 *        block
 492          IDONE = 0
 493 *        loop while( IDONE.LT.IM )
 494 *        generate the representation tree for the current block and
 495 *        compute the eigenvectors
 496    40    CONTINUE
 497          IF( IDONE.LT.IM ) THEN
 498 *           This is a crude protection against infinitely deep trees
 499             IF( NDEPTH.GT.M ) THEN
 500                INFO = -2
 501                RETURN
 502             ENDIF
 503 *           breadth first processing of the current level of the representation
 504 *           tree: OLDNCL = number of clusters on current level
 505             OLDNCL = NCLUS
 506 *           reset NCLUS to count the number of child clusters
 507             NCLUS = 0
 508 *
 509             PARITY = 1 - PARITY
 510             IF( PARITY.EQ.0 ) THEN
 511                OLDCLS = IINDC1
 512                NEWCLS = IINDC2
 513             ELSE
 514                OLDCLS = IINDC2
 515                NEWCLS = IINDC1
 516             END IF
 517 *           Process the clusters on the current level
 518             DO 150 I = 1, OLDNCL
 519                J = OLDCLS + 2*I
 520 *              OLDFST, OLDLST = first, last index of current cluster.
 521 *                               cluster indices start with 1 and are relative
 522 *                               to WBEGIN when accessing W, WGAP, WERR, Z
 523                OLDFST = IWORK( J-1 )
 524                OLDLST = IWORK( J )
 525                IF( NDEPTH.GT.0 ) THEN
 526 *                 Retrieve relatively robust representation (RRR) of cluster
 527 *                 that has been computed at the previous level
 528 *                 The RRR is stored in Z and overwritten once the eigenvectors
 529 *                 have been computed or when the cluster is refined
 530
 531                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
 532 *                    Get representation from location of the leftmost evalue
 533 *                    of the cluster
 534                      J = WBEGIN + OLDFST - 1
 535                   ELSE
 536                      IF(WBEGIN+OLDFST-1.LT.DOL) THEN
 537 *                       Get representation from the left end of Z array
 538                         J = DOL - 1
 539                      ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
 540 *                       Get representation from the right end of Z array
 541                         J = DOU
 542                      ELSE
 543                         J = WBEGIN + OLDFST - 1
 544                      ENDIF
 545                   ENDIF
 546                   CALL DCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 )
 547                   CALL DCOPY( IN-1, Z( IBEGIN, J+1 ), 1, L( IBEGIN ),
 548      $               1 )
 549                   SIGMA = Z( IEND, J+1 )
 550
 551 *                 Set the corresponding entries in Z to zero
 552                   CALL DLASET( 'Full', IN, 2, ZERO, ZERO,
 553      $                         Z( IBEGIN, J), LDZ )
 554                END IF
 555
 556 *              Compute DL and DLL of current RRR
 557                DO 50 J = IBEGIN, IEND-1
 558                   TMP = D( J )*L( J )
 559                   WORK( INDLD-1+J ) = TMP
 560                   WORK( INDLLD-1+J ) = TMP*L( J )
 561    50          CONTINUE
 562
 563                IF( NDEPTH.GT.0 ) THEN
 564 *                 P and Q are index of the first and last eigenvalue to compute
 565 *                 within the current block
 566                   P = INDEXW( WBEGIN-1+OLDFST )
 567                   Q = INDEXW( WBEGIN-1+OLDLST )
 568 *                 Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
 569 *                 through the Q-OFFSET elements of these arrays are to be used.
 570 *                  OFFSET = P-OLDFST
 571                   OFFSET = INDEXW( WBEGIN ) - 1
 572 *                 perform limited bisection (if necessary) to get approximate
 573 *                 eigenvalues to the precision needed.
 574                   CALL DLARRB( IN, D( IBEGIN ),
 575      $                         WORK(INDLLD+IBEGIN-1),
 576      $                         P, Q, RTOL1, RTOL2, OFFSET,
 577      $                         WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
 578      $                         WORK( INDWRK ), IWORK( IINDWK ),
 579      $                         PIVMIN, SPDIAM, IN, IINFO )
 580                   IF( IINFO.NE.0 ) THEN
 581                      INFO = -1
 582                      RETURN
 583                   ENDIF
 584 *                 We also recompute the extremal gaps. W holds all eigenvalues
 585 *                 of the unshifted matrix and must be used for computation
 586 *                 of WGAP, the entries of WORK might stem from RRRs with
 587 *                 different shifts. The gaps from WBEGIN-1+OLDFST to
 588 *                 WBEGIN-1+OLDLST are correctly computed in DLARRB.
 589 *                 However, we only allow the gaps to become greater since
 590 *                 this is what should happen when we decrease WERR
 591                   IF( OLDFST.GT.1) THEN
 592                      WGAP( WBEGIN+OLDFST-2 ) =
 593      $             MAX(WGAP(WBEGIN+OLDFST-2),
 594      $                 W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
 595      $                 - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
 596                   ENDIF
 597                   IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
 598                      WGAP( WBEGIN+OLDLST-1 ) =
 599      $               MAX(WGAP(WBEGIN+OLDLST-1),
 600      $                   W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
 601      $                   - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
 602                   ENDIF
 603 *                 Each time the eigenvalues in WORK get refined, we store
 604 *                 the newly found approximation with all shifts applied in W
 605                   DO 53 J=OLDFST,OLDLST
 606                      W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
 607  53               CONTINUE
 608                END IF
 609
 610 *              Process the current node.
 611                NEWFST = OLDFST
 612                DO 140 J = OLDFST, OLDLST
 613                   IF( J.EQ.OLDLST ) THEN
 614 *                    we are at the right end of the cluster, this is also the
 615 *                    boundary of the child cluster
 616                      NEWLST = J
 617                   ELSE IF ( WGAP( WBEGIN + J -1).GE.
 618      $                    MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
 619 *                    the right relative gap is big enough, the child cluster
 620 *                    (NEWFST,..,NEWLST) is well separated from the following
 621                      NEWLST = J
 622                    ELSE
 623 *                    inside a child cluster, the relative gap is not
 624 *                    big enough.
 625                      GOTO 140
 626                   END IF
 627
 628 *                 Compute size of child cluster found
 629                   NEWSIZ = NEWLST - NEWFST + 1
 630
 631 *                 NEWFTT is the place in Z where the new RRR or the computed
 632 *                 eigenvector is to be stored
 633                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
 634 *                    Store representation at location of the leftmost evalue
 635 *                    of the cluster
 636                      NEWFTT = WBEGIN + NEWFST - 1
 637                   ELSE
 638                      IF(WBEGIN+NEWFST-1.LT.DOL) THEN
 639 *                       Store representation at the left end of Z array
 640                         NEWFTT = DOL - 1
 641                      ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
 642 *                       Store representation at the right end of Z array
 643                         NEWFTT = DOU
 644                      ELSE
 645                         NEWFTT = WBEGIN + NEWFST - 1
 646                      ENDIF
 647                   ENDIF
 648
 649                   IF( NEWSIZ.GT.1) THEN
 650 *
 651 *                    Current child is not a singleton but a cluster.
 652 *                    Compute and store new representation of child.
 653 *
 654 *
 655 *                    Compute left and right cluster gap.
 656 *
 657 *                    LGAP and RGAP are not computed from WORK because
 658 *                    the eigenvalue approximations may stem from RRRs
 659 *                    different shifts. However, W hold all eigenvalues
 660 *                    of the unshifted matrix. Still, the entries in WGAP
 661 *                    have to be computed from WORK since the entries
 662 *                    in W might be of the same order so that gaps are not
 663 *                    exhibited correctly for very close eigenvalues.
 664                      IF( NEWFST.EQ.1 ) THEN
 665                         LGAP = MAX( ZERO,
 666      $                       W(WBEGIN)-WERR(WBEGIN) - VL )
 667                     ELSE
 668                         LGAP = WGAP( WBEGIN+NEWFST-2 )
 669                      ENDIF
 670                      RGAP = WGAP( WBEGIN+NEWLST-1 )
 671 *
 672 *                    Compute left- and rightmost eigenvalue of child
 673 *                    to high precision in order to shift as close
 674 *                    as possible and obtain as large relative gaps
 675 *                    as possible
 676 *
 677                      DO 55 K =1,2
 678                         IF(K.EQ.1) THEN
 679                            P = INDEXW( WBEGIN-1+NEWFST )
 680                         ELSE
 681                            P = INDEXW( WBEGIN-1+NEWLST )
 682                         ENDIF
 683                         OFFSET = INDEXW( WBEGIN ) - 1
 684                         CALL DLARRB( IN, D(IBEGIN),
 685      $                       WORK( INDLLD+IBEGIN-1 ),P,P,
 686      $                       RQTOL, RQTOL, OFFSET,
 687      $                       WORK(WBEGIN),WGAP(WBEGIN),
 688      $                       WERR(WBEGIN),WORK( INDWRK ),
 689      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
 690      $                       IN, IINFO )
 691  55                  CONTINUE
 692 *
 693                      IF((WBEGIN+NEWLST-1.LT.DOL).OR.
 694      $                  (WBEGIN+NEWFST-1.GT.DOU)) THEN
 695 *                       if the cluster contains no desired eigenvalues
 696 *                       skip the computation of that branch of the rep. tree
 697 *
 698 *                       We could skip before the refinement of the extremal
 699 *                       eigenvalues of the child, but then the representation
 700 *                       tree could be different from the one when nothing is
 701 *                       skipped. For this reason we skip at this place.
 702                         IDONE = IDONE + NEWLST - NEWFST + 1
 703                         GOTO 139
 704                      ENDIF
 705 *
 706 *                    Compute RRR of child cluster.
 707 *                    Note that the new RRR is stored in Z
 708 *
 709 *                    DLARRF needs LWORK = 2*N
 710                      CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ),
 711      $                         WORK(INDLD+IBEGIN-1),
 712      $                         NEWFST, NEWLST, WORK(WBEGIN),
 713      $                         WGAP(WBEGIN), WERR(WBEGIN),
 714      $                         SPDIAM, LGAP, RGAP, PIVMIN, TAU,
 715      $                         Z(IBEGIN, NEWFTT),Z(IBEGIN, NEWFTT+1),
 716      $                         WORK( INDWRK ), IINFO )
 717                      IF( IINFO.EQ.0 ) THEN
 718 *                       a new RRR for the cluster was found by DLARRF
 719 *                       update shift and store it
 720                         SSIGMA = SIGMA + TAU
 721                         Z( IEND, NEWFTT+1 ) = SSIGMA
 722 *                       WORK() are the midpoints and WERR() the semi-width
 723 *                       Note that the entries in W are unchanged.
 724                         DO 116 K = NEWFST, NEWLST
 725                            FUDGE =
 726      $                          THREE*EPS*ABS(WORK(WBEGIN+K-1))
 727                            WORK( WBEGIN + K - 1 ) =
 728      $                          WORK( WBEGIN + K - 1) - TAU
 729                            FUDGE = FUDGE +
 730      $                          FOUR*EPS*ABS(WORK(WBEGIN+K-1))
 731 *                          Fudge errors
 732                            WERR( WBEGIN + K - 1 ) =
 733      $                          WERR( WBEGIN + K - 1 ) + FUDGE
 734 *                          Gaps are not fudged. Provided that WERR is small
 735 *                          when eigenvalues are close, a zero gap indicates
 736 *                          that a new representation is needed for resolving
 737 *                          the cluster. A fudge could lead to a wrong decision
 738 *                          of judging eigenvalues 'separated' which in
 739 *                          reality are not. This could have a negative impact
 740 *                          on the orthogonality of the computed eigenvectors.
 741  116                    CONTINUE
 742
 743                         NCLUS = NCLUS + 1
 744                         K = NEWCLS + 2*NCLUS
 745                         IWORK( K-1 ) = NEWFST
 746                         IWORK( K ) = NEWLST
 747                      ELSE
 748                         INFO = -2
 749                         RETURN
 750                      ENDIF
 751                   ELSE
 752 *
 753 *                    Compute eigenvector of singleton
 754 *
 755                      ITER = 0
 756 *
 757                      TOL = FOUR * LOG(DBLE(IN)) * EPS
 758 *
 759                      K = NEWFST
 760                      WINDEX = WBEGIN + K - 1
 761                      WINDMN = MAX(WINDEX - 1,1)
 762                      WINDPL = MIN(WINDEX + 1,M)
 763                      LAMBDA = WORK( WINDEX )
 764                      DONE = DONE + 1
 765 *                    Check if eigenvector computation is to be skipped
 766                      IF((WINDEX.LT.DOL).OR.
 767      $                  (WINDEX.GT.DOU)) THEN
 768                         ESKIP = .TRUE.
 769                         GOTO 125
 770                      ELSE
 771                         ESKIP = .FALSE.
 772                      ENDIF
 773                      LEFT = WORK( WINDEX ) - WERR( WINDEX )
 774                      RIGHT = WORK( WINDEX ) + WERR( WINDEX )
 775                      INDEIG = INDEXW( WINDEX )
 776 *                    Note that since we compute the eigenpairs for a child,
 777 *                    all eigenvalue approximations are w.r.t the same shift.
 778 *                    In this case, the entries in WORK should be used for
 779 *                    computing the gaps since they exhibit even very small
 780 *                    differences in the eigenvalues, as opposed to the
 781 *                    entries in W which might "look" the same.
 782
 783                      IF( K .EQ. 1) THEN
 784 *                       In the case RANGE='I' and with not much initial
 785 *                       accuracy in LAMBDA and VL, the formula
 786 *                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
 787 *                       can lead to an overestimation of the left gap and
 788 *                       thus to inadequately early RQI 'convergence'.
 789 *                       Prevent this by forcing a small left gap.
 790                         LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
 791                      ELSE
 792                         LGAP = WGAP(WINDMN)
 793                      ENDIF
 794                      IF( K .EQ. IM) THEN
 795 *                       In the case RANGE='I' and with not much initial
 796 *                       accuracy in LAMBDA and VU, the formula
 797 *                       can lead to an overestimation of the right gap and
 798 *                       thus to inadequately early RQI 'convergence'.
 799 *                       Prevent this by forcing a small right gap.
 800                         RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
 801                      ELSE
 802                         RGAP = WGAP(WINDEX)
 803                      ENDIF
 804                      GAP = MIN( LGAP, RGAP )
 805                      IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN
 806 *                       The eigenvector support can become wrong
 807 *                       because significant entries could be cut off due to a
 808 *                       large GAPTOL parameter in LAR1V. Prevent this.
 809                         GAPTOL = ZERO
 810                      ELSE
 811                         GAPTOL = GAP * EPS
 812                      ENDIF
 813                      ISUPMN = IN
 814                      ISUPMX = 1
 815 *                    Update WGAP so that it holds the minimum gap
 816 *                    to the left or the right. This is crucial in the
 817 *                    case where bisection is used to ensure that the
 818 *                    eigenvalue is refined up to the required precision.
 819 *                    The correct value is restored afterwards.
 820                      SAVGAP = WGAP(WINDEX)
 821                      WGAP(WINDEX) = GAP
 822 *                    We want to use the Rayleigh Quotient Correction
 823 *                    as often as possible since it converges quadratically
 824 *                    when we are close enough to the desired eigenvalue.
 825 *                    However, the Rayleigh Quotient can have the wrong sign
 826 *                    and lead us away from the desired eigenvalue. In this
 827 *                    case, the best we can do is to use bisection.
 828                      USEDBS = .FALSE.
 829                      USEDRQ = .FALSE.
 830 *                    Bisection is initially turned off unless it is forced
 831                      NEEDBS =  .NOT.TRYRQC
 832  120                 CONTINUE
 833 *                    Check if bisection should be used to refine eigenvalue
 834                      IF(NEEDBS) THEN
 835 *                       Take the bisection as new iterate
 836                         USEDBS = .TRUE.
 837                         ITMP1 = IWORK( IINDR+WINDEX )
 838                         OFFSET = INDEXW( WBEGIN ) - 1
 839                         CALL DLARRB( IN, D(IBEGIN),
 840      $                       WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
 841      $                       ZERO, TWO*EPS, OFFSET,
 842      $                       WORK(WBEGIN),WGAP(WBEGIN),
 843      $                       WERR(WBEGIN),WORK( INDWRK ),
 844      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
 845      $                       ITMP1, IINFO )
 846                         IF( IINFO.NE.0 ) THEN
 847                            INFO = -3
 848                            RETURN
 849                         ENDIF
 850                         LAMBDA = WORK( WINDEX )
 851 *                       Reset twist index from inaccurate LAMBDA to
 852 *                       force computation of true MINGMA
 853                         IWORK( IINDR+WINDEX ) = 0
 854                      ENDIF
 855 *                    Given LAMBDA, compute the eigenvector.
 856                      CALL DLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ),
 857      $                    L( IBEGIN ), WORK(INDLD+IBEGIN-1),
 858      $                    WORK(INDLLD+IBEGIN-1),
 859      $                    PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
 860      $                    .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
 861      $                    IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
 862      $                    NRMINV, RESID, RQCORR, WORK( INDWRK ) )
 863                      IF(ITER .EQ. 0) THEN
 864                         BSTRES = RESID
 865                         BSTW = LAMBDA
 866                      ELSEIF(RESID.LT.BSTRES) THEN
 867                         BSTRES = RESID
 868                         BSTW = LAMBDA
 869                      ENDIF
 870                      ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
 871                      ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
 872                      ITER = ITER + 1
 873
 874 *                    sin alpha <= |resid|/gap
 875 *                    Note that both the residual and the gap are
 876 *                    proportional to the matrix, so ||T|| doesn't play
 877 *                    a role in the quotient
 878
 879 *
 880 *                    Convergence test for Rayleigh-Quotient iteration
 881 *                    (omitted when Bisection has been used)
 882 *
 883                      IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
 884      $                    RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
 885      $                    THEN
 886 *                       We need to check that the RQCORR update doesn't
 887 *                       move the eigenvalue away from the desired one and
 888 *                       towards a neighbor. -> protection with bisection
 889                         IF(INDEIG.LE.NEGCNT) THEN
 890 *                          The wanted eigenvalue lies to the left
 891                            SGNDEF = -ONE
 892                         ELSE
 893 *                          The wanted eigenvalue lies to the right
 894                            SGNDEF = ONE
 895                         ENDIF
 896 *                       We only use the RQCORR if it improves the
 897 *                       the iterate reasonably.
 898                         IF( ( RQCORR*SGNDEF.GE.ZERO )
 899      $                       .AND.( LAMBDA + RQCORR.LE. RIGHT)
 900      $                       .AND.( LAMBDA + RQCORR.GE. LEFT)
 901      $                       ) THEN
 902                            USEDRQ = .TRUE.
 903 *                          Store new midpoint of bisection interval in WORK
 904                            IF(SGNDEF.EQ.ONE) THEN
 905 *                             The current LAMBDA is on the left of the true
 906 *                             eigenvalue
 907                               LEFT = LAMBDA
 908 *                             We prefer to assume that the error estimate
 909 *                             is correct. We could make the interval not
 910 *                             as a bracket but to be modified if the RQCORR
 911 *                             chooses to. In this case, the RIGHT side should
 912 *                             be modified as follows:
 913 *                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
 914                            ELSE
 915 *                             The current LAMBDA is on the right of the true
 916 *                             eigenvalue
 917                               RIGHT = LAMBDA
 918 *                             See comment about assuming the error estimate is
 919 *                             correct above.
 920 *                              LEFT = MIN(LEFT, LAMBDA + RQCORR)
 921                            ENDIF
 922                            WORK( WINDEX ) =
 923      $                       HALF * (RIGHT + LEFT)
 924 *                          Take RQCORR since it has the correct sign and
 925 *                          improves the iterate reasonably
 926                            LAMBDA = LAMBDA + RQCORR
 927 *                          Update width of error interval
 928                            WERR( WINDEX ) =
 929      $                             HALF * (RIGHT-LEFT)
 930                         ELSE
 931                            NEEDBS = .TRUE.
 932                         ENDIF
 933                         IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
 934 *                             The eigenvalue is computed to bisection accuracy
 935 *                             compute eigenvector and stop
 936                            USEDBS = .TRUE.
 937                            GOTO 120
 938                         ELSEIF( ITER.LT.MAXITR ) THEN
 939                            GOTO 120
 940                         ELSEIF( ITER.EQ.MAXITR ) THEN
 941                            NEEDBS = .TRUE.
 942                            GOTO 120
 943                         ELSE
 944                            INFO = 5
 945                            RETURN
 946                         END IF
 947                      ELSE
 948                         STP2II = .FALSE.
 949         IF(USEDRQ .AND. USEDBS .AND.
 950      $                     BSTRES.LE.RESID) THEN
 951                            LAMBDA = BSTW
 952                            STP2II = .TRUE.
 953                         ENDIF
 954                         IF (STP2II) THEN
 955 *                          improve error angle by second step
 956                            CALL DLAR1V( IN, 1, IN, LAMBDA,
 957      $                          D( IBEGIN ), L( IBEGIN ),
 958      $                          WORK(INDLD+IBEGIN-1),
 959      $                          WORK(INDLLD+IBEGIN-1),
 960      $                          PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
 961      $                          .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
 962      $                          IWORK( IINDR+WINDEX ),
 963      $                          ISUPPZ( 2*WINDEX-1 ),
 964      $                          NRMINV, RESID, RQCORR, WORK( INDWRK ) )
 965                         ENDIF
 966                         WORK( WINDEX ) = LAMBDA
 967                      END IF
 968 *
 969 *                    Compute FP-vector support w.r.t. whole matrix
 970 *
 971                      ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
 972                      ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
 973                      ZFROM = ISUPPZ( 2*WINDEX-1 )
 974                      ZTO = ISUPPZ( 2*WINDEX )
 975                      ISUPMN = ISUPMN + OLDIEN
 976                      ISUPMX = ISUPMX + OLDIEN
 977 *                    Ensure vector is ok if support in the RQI has changed
 978                      IF(ISUPMN.LT.ZFROM) THEN
 979                         DO 122 II = ISUPMN,ZFROM-1
 980                            Z( II, WINDEX ) = ZERO
 981  122                    CONTINUE
 982                      ENDIF
 983                      IF(ISUPMX.GT.ZTO) THEN
 984                         DO 123 II = ZTO+1,ISUPMX
 985                            Z( II, WINDEX ) = ZERO
 986  123                    CONTINUE
 987                      ENDIF
 988                      CALL DSCAL( ZTO-ZFROM+1, NRMINV,
 989      $                       Z( ZFROM, WINDEX ), 1 )
 990  125                 CONTINUE
 991 *                    Update W
 992                      W( WINDEX ) = LAMBDA+SIGMA
 993 *                    Recompute the gaps on the left and right
 994 *                    But only allow them to become larger and not
 995 *                    smaller (which can only happen through "bad"
 996 *                    cancellation and doesn't reflect the theory
 997 *                    where the initial gaps are underestimated due
 998 *                    to WERR being too crude.)
 999                      IF(.NOT.ESKIP) THEN
1000                         IF( K.GT.1) THEN
1001                            WGAP( WINDMN ) = MAX( WGAP(WINDMN),
1002      $                          W(WINDEX)-WERR(WINDEX)
1003      $                          - W(WINDMN)-WERR(WINDMN) )
1004                         ENDIF
1005                         IF( WINDEX.LT.WEND ) THEN
1006                            WGAP( WINDEX ) = MAX( SAVGAP,
1007      $                          W( WINDPL )-WERR( WINDPL )
1008      $                          - W( WINDEX )-WERR( WINDEX) )
1009                         ENDIF
1010                      ENDIF
1011                      IDONE = IDONE + 1
1012                   ENDIF
1013 *                 here ends the code for the current child
1014 *
1015  139              CONTINUE
1016 *                 Proceed to any remaining child nodes
1017                   NEWFST = J + 1
1018  140           CONTINUE
1019  150        CONTINUE
1020             NDEPTH = NDEPTH + 1
1021             GO TO 40
1022          END IF
1023          IBEGIN = IEND + 1
1024          WBEGIN = WEND + 1
1025  170  CONTINUE
1026 *
1027
1028       RETURN
1029 *
1030 *     End of DLARRV
1031 *
1032       END