SRC/zlarrv.f

   1 *> \brief \b ZLARRV 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 ZLARRV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarrv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarrv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarrv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLARRV( 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 *       COMPLEX*16        Z( LDZ, * )
  37 *       ..
  38 *
  39 *
  40 *> \par Purpose:
  41 *  =============
  42 *>
  43 *> \verbatim
  44 *>
  45 *> ZLARRV 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 COMPLEX*16 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 ZLARRV.
 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 complex16OTHERauxiliary
 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 ZLARRV( 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       COMPLEX*16        Z( LDZ, * )
 304 *     ..
 305 *
 306 *  =====================================================================
 307 *
 308 *     .. Parameters ..
 309       INTEGER            MAXITR
 310       PARAMETER          ( MAXITR = 10 )
 311       COMPLEX*16         CZERO
 312       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
 313       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, HALF
 314       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
 315      $                     TWO = 2.0D0, THREE = 3.0D0,
 316      $                     FOUR = 4.0D0, HALF = 0.5D0)
 317 *     ..
 318 *     .. Local Scalars ..
 319       LOGICAL            ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
 320       INTEGER            DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
 321      $                   IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
 322      $                   INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
 323      $                   ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
 324      $                   NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
 325      $                   NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
 326      $                   OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
 327      $                   WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
 328      $                   ZUSEDW
 329       INTEGER            INDIN1, INDIN2
 330       DOUBLE PRECISION   BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
 331      $                   LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
 332      $                   RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
 333      $                   SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
 334 *     ..
 335 *     .. External Functions ..
 336       DOUBLE PRECISION   DLAMCH
 337       EXTERNAL           DLAMCH
 338 *     ..
 339 *     .. External Subroutines ..
 340       EXTERNAL           DCOPY, DLARRB, DLARRF, ZDSCAL, ZLAR1V,
 341      $                   ZLASET
 342 *     ..
 343 *     .. Intrinsic Functions ..
 344       INTRINSIC ABS, DBLE, MAX, MIN
 345       INTRINSIC DCMPLX
 346 *     ..
 347 *     .. Executable Statements ..
 348 *     ..
 349
 350       INFO = 0
 351 *     The first N entries of WORK are reserved for the eigenvalues
 352       INDLD = N+1
 353       INDLLD= 2*N+1
 354       INDIN1 = 3*N + 1
 355       INDIN2 = 4*N + 1
 356       INDWRK = 5*N + 1
 357       MINWSIZE = 12 * N
 358
 359       DO 5 I= 1,MINWSIZE
 360          WORK( I ) = ZERO
 361  5    CONTINUE
 362
 363 *     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
 364 *     factorization used to compute the FP vector
 365       IINDR = 0
 366 *     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
 367 *     layer and the one above.
 368       IINDC1 = N
 369       IINDC2 = 2*N
 370       IINDWK = 3*N + 1
 371
 372       MINIWSIZE = 7 * N
 373       DO 10 I= 1,MINIWSIZE
 374          IWORK( I ) = 0
 375  10   CONTINUE
 376
 377       ZUSEDL = 1
 378       IF(DOL.GT.1) THEN
 379 *        Set lower bound for use of Z
 380          ZUSEDL = DOL-1
 381       ENDIF
 382       ZUSEDU = M
 383       IF(DOU.LT.M) THEN
 384 *        Set lower bound for use of Z
 385          ZUSEDU = DOU+1
 386       ENDIF
 387 *     The width of the part of Z that is used
 388       ZUSEDW = ZUSEDU - ZUSEDL + 1
 389
 390
 391       CALL ZLASET( 'Full', N, ZUSEDW, CZERO, CZERO,
 392      $                    Z(1,ZUSEDL), LDZ )
 393
 394       EPS = DLAMCH( 'Precision' )
 395       RQTOL = TWO * EPS
 396 *
 397 *     Set expert flags for standard code.
 398       TRYRQC = .TRUE.
 399
 400       IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
 401       ELSE
 402 *        Only selected eigenpairs are computed. Since the other evalues
 403 *        are not refined by RQ iteration, bisection has to compute to full
 404 *        accuracy.
 405          RTOL1 = FOUR * EPS
 406          RTOL2 = FOUR * EPS
 407       ENDIF
 408
 409 *     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
 410 *     desired eigenvalues. The support of the nonzero eigenvector
 411 *     entries is contained in the interval IBEGIN:IEND.
 412 *     Remark that if k eigenpairs are desired, then the eigenvectors
 413 *     are stored in k contiguous columns of Z.
 414
 415 *     DONE is the number of eigenvectors already computed
 416       DONE = 0
 417       IBEGIN = 1
 418       WBEGIN = 1
 419       DO 170 JBLK = 1, IBLOCK( M )
 420          IEND = ISPLIT( JBLK )
 421          SIGMA = L( IEND )
 422 *        Find the eigenvectors of the submatrix indexed IBEGIN
 423 *        through IEND.
 424          WEND = WBEGIN - 1
 425  15      CONTINUE
 426          IF( WEND.LT.M ) THEN
 427             IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
 428                WEND = WEND + 1
 429                GO TO 15
 430             END IF
 431          END IF
 432          IF( WEND.LT.WBEGIN ) THEN
 433             IBEGIN = IEND + 1
 434             GO TO 170
 435          ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
 436             IBEGIN = IEND + 1
 437             WBEGIN = WEND + 1
 438             GO TO 170
 439          END IF
 440
 441 *        Find local spectral diameter of the block
 442          GL = GERS( 2*IBEGIN-1 )
 443          GU = GERS( 2*IBEGIN )
 444          DO 20 I = IBEGIN+1 , IEND
 445             GL = MIN( GERS( 2*I-1 ), GL )
 446             GU = MAX( GERS( 2*I ), GU )
 447  20      CONTINUE
 448          SPDIAM = GU - GL
 449
 450 *        OLDIEN is the last index of the previous block
 451          OLDIEN = IBEGIN - 1
 452 *        Calculate the size of the current block
 453          IN = IEND - IBEGIN + 1
 454 *        The number of eigenvalues in the current block
 455          IM = WEND - WBEGIN + 1
 456
 457 *        This is for a 1x1 block
 458          IF( IBEGIN.EQ.IEND ) THEN
 459             DONE = DONE+1
 460             Z( IBEGIN, WBEGIN ) = DCMPLX( ONE, ZERO )
 461             ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
 462             ISUPPZ( 2*WBEGIN ) = IBEGIN
 463             W( WBEGIN ) = W( WBEGIN ) + SIGMA
 464             WORK( WBEGIN ) = W( WBEGIN )
 465             IBEGIN = IEND + 1
 466             WBEGIN = WBEGIN + 1
 467             GO TO 170
 468          END IF
 469
 470 *        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
 471 *        Note that these can be approximations, in this case, the corresp.
 472 *        entries of WERR give the size of the uncertainty interval.
 473 *        The eigenvalue approximations will be refined when necessary as
 474 *        high relative accuracy is required for the computation of the
 475 *        corresponding eigenvectors.
 476          CALL DCOPY( IM, W( WBEGIN ), 1,
 477      $                   WORK( WBEGIN ), 1 )
 478
 479 *        We store in W the eigenvalue approximations w.r.t. the original
 480 *        matrix T.
 481          DO 30 I=1,IM
 482             W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
 483  30      CONTINUE
 484
 485
 486 *        NDEPTH is the current depth of the representation tree
 487          NDEPTH = 0
 488 *        PARITY is either 1 or 0
 489          PARITY = 1
 490 *        NCLUS is the number of clusters for the next level of the
 491 *        representation tree, we start with NCLUS = 1 for the root
 492          NCLUS = 1
 493          IWORK( IINDC1+1 ) = 1
 494          IWORK( IINDC1+2 ) = IM
 495
 496 *        IDONE is the number of eigenvectors already computed in the current
 497 *        block
 498          IDONE = 0
 499 *        loop while( IDONE.LT.IM )
 500 *        generate the representation tree for the current block and
 501 *        compute the eigenvectors
 502    40    CONTINUE
 503          IF( IDONE.LT.IM ) THEN
 504 *           This is a crude protection against infinitely deep trees
 505             IF( NDEPTH.GT.M ) THEN
 506                INFO = -2
 507                RETURN
 508             ENDIF
 509 *           breadth first processing of the current level of the representation
 510 *           tree: OLDNCL = number of clusters on current level
 511             OLDNCL = NCLUS
 512 *           reset NCLUS to count the number of child clusters
 513             NCLUS = 0
 514 *
 515             PARITY = 1 - PARITY
 516             IF( PARITY.EQ.0 ) THEN
 517                OLDCLS = IINDC1
 518                NEWCLS = IINDC2
 519             ELSE
 520                OLDCLS = IINDC2
 521                NEWCLS = IINDC1
 522             END IF
 523 *           Process the clusters on the current level
 524             DO 150 I = 1, OLDNCL
 525                J = OLDCLS + 2*I
 526 *              OLDFST, OLDLST = first, last index of current cluster.
 527 *                               cluster indices start with 1 and are relative
 528 *                               to WBEGIN when accessing W, WGAP, WERR, Z
 529                OLDFST = IWORK( J-1 )
 530                OLDLST = IWORK( J )
 531                IF( NDEPTH.GT.0 ) THEN
 532 *                 Retrieve relatively robust representation (RRR) of cluster
 533 *                 that has been computed at the previous level
 534 *                 The RRR is stored in Z and overwritten once the eigenvectors
 535 *                 have been computed or when the cluster is refined
 536
 537                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
 538 *                    Get representation from location of the leftmost evalue
 539 *                    of the cluster
 540                      J = WBEGIN + OLDFST - 1
 541                   ELSE
 542                      IF(WBEGIN+OLDFST-1.LT.DOL) THEN
 543 *                       Get representation from the left end of Z array
 544                         J = DOL - 1
 545                      ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
 546 *                       Get representation from the right end of Z array
 547                         J = DOU
 548                      ELSE
 549                         J = WBEGIN + OLDFST - 1
 550                      ENDIF
 551                   ENDIF
 552                   DO 45 K = 1, IN - 1
 553                      D( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1,
 554      $                                 J ) )
 555                      L( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1,
 556      $                                 J+1 ) )
 557    45             CONTINUE
 558                   D( IEND ) = DBLE( Z( IEND, J ) )
 559                   SIGMA = DBLE( Z( IEND, J+1 ) )
 560
 561 *                 Set the corresponding entries in Z to zero
 562                   CALL ZLASET( 'Full', IN, 2, CZERO, CZERO,
 563      $                         Z( IBEGIN, J), LDZ )
 564                END IF
 565
 566 *              Compute DL and DLL of current RRR
 567                DO 50 J = IBEGIN, IEND-1
 568                   TMP = D( J )*L( J )
 569                   WORK( INDLD-1+J ) = TMP
 570                   WORK( INDLLD-1+J ) = TMP*L( J )
 571    50          CONTINUE
 572
 573                IF( NDEPTH.GT.0 ) THEN
 574 *                 P and Q are index of the first and last eigenvalue to compute
 575 *                 within the current block
 576                   P = INDEXW( WBEGIN-1+OLDFST )
 577                   Q = INDEXW( WBEGIN-1+OLDLST )
 578 *                 Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
 579 *                 through the Q-OFFSET elements of these arrays are to be used.
 580 *                  OFFSET = P-OLDFST
 581                   OFFSET = INDEXW( WBEGIN ) - 1
 582 *                 perform limited bisection (if necessary) to get approximate
 583 *                 eigenvalues to the precision needed.
 584                   CALL DLARRB( IN, D( IBEGIN ),
 585      $                         WORK(INDLLD+IBEGIN-1),
 586      $                         P, Q, RTOL1, RTOL2, OFFSET,
 587      $                         WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
 588      $                         WORK( INDWRK ), IWORK( IINDWK ),
 589      $                         PIVMIN, SPDIAM, IN, IINFO )
 590                   IF( IINFO.NE.0 ) THEN
 591                      INFO = -1
 592                      RETURN
 593                   ENDIF
 594 *                 We also recompute the extremal gaps. W holds all eigenvalues
 595 *                 of the unshifted matrix and must be used for computation
 596 *                 of WGAP, the entries of WORK might stem from RRRs with
 597 *                 different shifts. The gaps from WBEGIN-1+OLDFST to
 598 *                 WBEGIN-1+OLDLST are correctly computed in DLARRB.
 599 *                 However, we only allow the gaps to become greater since
 600 *                 this is what should happen when we decrease WERR
 601                   IF( OLDFST.GT.1) THEN
 602                      WGAP( WBEGIN+OLDFST-2 ) =
 603      $             MAX(WGAP(WBEGIN+OLDFST-2),
 604      $                 W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
 605      $                 - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
 606                   ENDIF
 607                   IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
 608                      WGAP( WBEGIN+OLDLST-1 ) =
 609      $               MAX(WGAP(WBEGIN+OLDLST-1),
 610      $                   W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
 611      $                   - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
 612                   ENDIF
 613 *                 Each time the eigenvalues in WORK get refined, we store
 614 *                 the newly found approximation with all shifts applied in W
 615                   DO 53 J=OLDFST,OLDLST
 616                      W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
 617  53               CONTINUE
 618                END IF
 619
 620 *              Process the current node.
 621                NEWFST = OLDFST
 622                DO 140 J = OLDFST, OLDLST
 623                   IF( J.EQ.OLDLST ) THEN
 624 *                    we are at the right end of the cluster, this is also the
 625 *                    boundary of the child cluster
 626                      NEWLST = J
 627                   ELSE IF ( WGAP( WBEGIN + J -1).GE.
 628      $                    MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
 629 *                    the right relative gap is big enough, the child cluster
 630 *                    (NEWFST,..,NEWLST) is well separated from the following
 631                      NEWLST = J
 632                    ELSE
 633 *                    inside a child cluster, the relative gap is not
 634 *                    big enough.
 635                      GOTO 140
 636                   END IF
 637
 638 *                 Compute size of child cluster found
 639                   NEWSIZ = NEWLST - NEWFST + 1
 640
 641 *                 NEWFTT is the place in Z where the new RRR or the computed
 642 *                 eigenvector is to be stored
 643                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
 644 *                    Store representation at location of the leftmost evalue
 645 *                    of the cluster
 646                      NEWFTT = WBEGIN + NEWFST - 1
 647                   ELSE
 648                      IF(WBEGIN+NEWFST-1.LT.DOL) THEN
 649 *                       Store representation at the left end of Z array
 650                         NEWFTT = DOL - 1
 651                      ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
 652 *                       Store representation at the right end of Z array
 653                         NEWFTT = DOU
 654                      ELSE
 655                         NEWFTT = WBEGIN + NEWFST - 1
 656                      ENDIF
 657                   ENDIF
 658
 659                   IF( NEWSIZ.GT.1) THEN
 660 *
 661 *                    Current child is not a singleton but a cluster.
 662 *                    Compute and store new representation of child.
 663 *
 664 *
 665 *                    Compute left and right cluster gap.
 666 *
 667 *                    LGAP and RGAP are not computed from WORK because
 668 *                    the eigenvalue approximations may stem from RRRs
 669 *                    different shifts. However, W hold all eigenvalues
 670 *                    of the unshifted matrix. Still, the entries in WGAP
 671 *                    have to be computed from WORK since the entries
 672 *                    in W might be of the same order so that gaps are not
 673 *                    exhibited correctly for very close eigenvalues.
 674                      IF( NEWFST.EQ.1 ) THEN
 675                         LGAP = MAX( ZERO,
 676      $                       W(WBEGIN)-WERR(WBEGIN) - VL )
 677                     ELSE
 678                         LGAP = WGAP( WBEGIN+NEWFST-2 )
 679                      ENDIF
 680                      RGAP = WGAP( WBEGIN+NEWLST-1 )
 681 *
 682 *                    Compute left- and rightmost eigenvalue of child
 683 *                    to high precision in order to shift as close
 684 *                    as possible and obtain as large relative gaps
 685 *                    as possible
 686 *
 687                      DO 55 K =1,2
 688                         IF(K.EQ.1) THEN
 689                            P = INDEXW( WBEGIN-1+NEWFST )
 690                         ELSE
 691                            P = INDEXW( WBEGIN-1+NEWLST )
 692                         ENDIF
 693                         OFFSET = INDEXW( WBEGIN ) - 1
 694                         CALL DLARRB( IN, D(IBEGIN),
 695      $                       WORK( INDLLD+IBEGIN-1 ),P,P,
 696      $                       RQTOL, RQTOL, OFFSET,
 697      $                       WORK(WBEGIN),WGAP(WBEGIN),
 698      $                       WERR(WBEGIN),WORK( INDWRK ),
 699      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
 700      $                       IN, IINFO )
 701  55                  CONTINUE
 702 *
 703                      IF((WBEGIN+NEWLST-1.LT.DOL).OR.
 704      $                  (WBEGIN+NEWFST-1.GT.DOU)) THEN
 705 *                       if the cluster contains no desired eigenvalues
 706 *                       skip the computation of that branch of the rep. tree
 707 *
 708 *                       We could skip before the refinement of the extremal
 709 *                       eigenvalues of the child, but then the representation
 710 *                       tree could be different from the one when nothing is
 711 *                       skipped. For this reason we skip at this place.
 712                         IDONE = IDONE + NEWLST - NEWFST + 1
 713                         GOTO 139
 714                      ENDIF
 715 *
 716 *                    Compute RRR of child cluster.
 717 *                    Note that the new RRR is stored in Z
 718 *
 719 *                    DLARRF needs LWORK = 2*N
 720                      CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ),
 721      $                         WORK(INDLD+IBEGIN-1),
 722      $                         NEWFST, NEWLST, WORK(WBEGIN),
 723      $                         WGAP(WBEGIN), WERR(WBEGIN),
 724      $                         SPDIAM, LGAP, RGAP, PIVMIN, TAU,
 725      $                         WORK( INDIN1 ), WORK( INDIN2 ),
 726      $                         WORK( INDWRK ), IINFO )
 727 *                    In the complex case, DLARRF cannot write
 728 *                    the new RRR directly into Z and needs an intermediate
 729 *                    workspace
 730                      DO 56 K = 1, IN-1
 731                         Z( IBEGIN+K-1, NEWFTT ) =
 732      $                     DCMPLX( WORK( INDIN1+K-1 ), ZERO )
 733                         Z( IBEGIN+K-1, NEWFTT+1 ) =
 734      $                     DCMPLX( WORK( INDIN2+K-1 ), ZERO )
 735    56                CONTINUE
 736                      Z( IEND, NEWFTT ) =
 737      $                  DCMPLX( WORK( INDIN1+IN-1 ), ZERO )
 738                      IF( IINFO.EQ.0 ) THEN
 739 *                       a new RRR for the cluster was found by DLARRF
 740 *                       update shift and store it
 741                         SSIGMA = SIGMA + TAU
 742                         Z( IEND, NEWFTT+1 ) = DCMPLX( SSIGMA, ZERO )
 743 *                       WORK() are the midpoints and WERR() the semi-width
 744 *                       Note that the entries in W are unchanged.
 745                         DO 116 K = NEWFST, NEWLST
 746                            FUDGE =
 747      $                          THREE*EPS*ABS(WORK(WBEGIN+K-1))
 748                            WORK( WBEGIN + K - 1 ) =
 749      $                          WORK( WBEGIN + K - 1) - TAU
 750                            FUDGE = FUDGE +
 751      $                          FOUR*EPS*ABS(WORK(WBEGIN+K-1))
 752 *                          Fudge errors
 753                            WERR( WBEGIN + K - 1 ) =
 754      $                          WERR( WBEGIN + K - 1 ) + FUDGE
 755 *                          Gaps are not fudged. Provided that WERR is small
 756 *                          when eigenvalues are close, a zero gap indicates
 757 *                          that a new representation is needed for resolving
 758 *                          the cluster. A fudge could lead to a wrong decision
 759 *                          of judging eigenvalues 'separated' which in
 760 *                          reality are not. This could have a negative impact
 761 *                          on the orthogonality of the computed eigenvectors.
 762  116                    CONTINUE
 763
 764                         NCLUS = NCLUS + 1
 765                         K = NEWCLS + 2*NCLUS
 766                         IWORK( K-1 ) = NEWFST
 767                         IWORK( K ) = NEWLST
 768                      ELSE
 769                         INFO = -2
 770                         RETURN
 771                      ENDIF
 772                   ELSE
 773 *
 774 *                    Compute eigenvector of singleton
 775 *
 776                      ITER = 0
 777 *
 778                      TOL = FOUR * LOG(DBLE(IN)) * EPS
 779 *
 780                      K = NEWFST
 781                      WINDEX = WBEGIN + K - 1
 782                      WINDMN = MAX(WINDEX - 1,1)
 783                      WINDPL = MIN(WINDEX + 1,M)
 784                      LAMBDA = WORK( WINDEX )
 785                      DONE = DONE + 1
 786 *                    Check if eigenvector computation is to be skipped
 787                      IF((WINDEX.LT.DOL).OR.
 788      $                  (WINDEX.GT.DOU)) THEN
 789                         ESKIP = .TRUE.
 790                         GOTO 125
 791                      ELSE
 792                         ESKIP = .FALSE.
 793                      ENDIF
 794                      LEFT = WORK( WINDEX ) - WERR( WINDEX )
 795                      RIGHT = WORK( WINDEX ) + WERR( WINDEX )
 796                      INDEIG = INDEXW( WINDEX )
 797 *                    Note that since we compute the eigenpairs for a child,
 798 *                    all eigenvalue approximations are w.r.t the same shift.
 799 *                    In this case, the entries in WORK should be used for
 800 *                    computing the gaps since they exhibit even very small
 801 *                    differences in the eigenvalues, as opposed to the
 802 *                    entries in W which might "look" the same.
 803
 804                      IF( K .EQ. 1) THEN
 805 *                       In the case RANGE='I' and with not much initial
 806 *                       accuracy in LAMBDA and VL, the formula
 807 *                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
 808 *                       can lead to an overestimation of the left gap and
 809 *                       thus to inadequately early RQI 'convergence'.
 810 *                       Prevent this by forcing a small left gap.
 811                         LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
 812                      ELSE
 813                         LGAP = WGAP(WINDMN)
 814                      ENDIF
 815                      IF( K .EQ. IM) THEN
 816 *                       In the case RANGE='I' and with not much initial
 817 *                       accuracy in LAMBDA and VU, the formula
 818 *                       can lead to an overestimation of the right gap and
 819 *                       thus to inadequately early RQI 'convergence'.
 820 *                       Prevent this by forcing a small right gap.
 821                         RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
 822                      ELSE
 823                         RGAP = WGAP(WINDEX)
 824                      ENDIF
 825                      GAP = MIN( LGAP, RGAP )
 826                      IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN
 827 *                       The eigenvector support can become wrong
 828 *                       because significant entries could be cut off due to a
 829 *                       large GAPTOL parameter in LAR1V. Prevent this.
 830                         GAPTOL = ZERO
 831                      ELSE
 832                         GAPTOL = GAP * EPS
 833                      ENDIF
 834                      ISUPMN = IN
 835                      ISUPMX = 1
 836 *                    Update WGAP so that it holds the minimum gap
 837 *                    to the left or the right. This is crucial in the
 838 *                    case where bisection is used to ensure that the
 839 *                    eigenvalue is refined up to the required precision.
 840 *                    The correct value is restored afterwards.
 841                      SAVGAP = WGAP(WINDEX)
 842                      WGAP(WINDEX) = GAP
 843 *                    We want to use the Rayleigh Quotient Correction
 844 *                    as often as possible since it converges quadratically
 845 *                    when we are close enough to the desired eigenvalue.
 846 *                    However, the Rayleigh Quotient can have the wrong sign
 847 *                    and lead us away from the desired eigenvalue. In this
 848 *                    case, the best we can do is to use bisection.
 849                      USEDBS = .FALSE.
 850                      USEDRQ = .FALSE.
 851 *                    Bisection is initially turned off unless it is forced
 852                      NEEDBS =  .NOT.TRYRQC
 853  120                 CONTINUE
 854 *                    Check if bisection should be used to refine eigenvalue
 855                      IF(NEEDBS) THEN
 856 *                       Take the bisection as new iterate
 857                         USEDBS = .TRUE.
 858                         ITMP1 = IWORK( IINDR+WINDEX )
 859                         OFFSET = INDEXW( WBEGIN ) - 1
 860                         CALL DLARRB( IN, D(IBEGIN),
 861      $                       WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
 862      $                       ZERO, TWO*EPS, OFFSET,
 863      $                       WORK(WBEGIN),WGAP(WBEGIN),
 864      $                       WERR(WBEGIN),WORK( INDWRK ),
 865      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
 866      $                       ITMP1, IINFO )
 867                         IF( IINFO.NE.0 ) THEN
 868                            INFO = -3
 869                            RETURN
 870                         ENDIF
 871                         LAMBDA = WORK( WINDEX )
 872 *                       Reset twist index from inaccurate LAMBDA to
 873 *                       force computation of true MINGMA
 874                         IWORK( IINDR+WINDEX ) = 0
 875                      ENDIF
 876 *                    Given LAMBDA, compute the eigenvector.
 877                      CALL ZLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ),
 878      $                    L( IBEGIN ), WORK(INDLD+IBEGIN-1),
 879      $                    WORK(INDLLD+IBEGIN-1),
 880      $                    PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
 881      $                    .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
 882      $                    IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
 883      $                    NRMINV, RESID, RQCORR, WORK( INDWRK ) )
 884                      IF(ITER .EQ. 0) THEN
 885                         BSTRES = RESID
 886                         BSTW = LAMBDA
 887                      ELSEIF(RESID.LT.BSTRES) THEN
 888                         BSTRES = RESID
 889                         BSTW = LAMBDA
 890                      ENDIF
 891                      ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
 892                      ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
 893                      ITER = ITER + 1
 894
 895 *                    sin alpha <= |resid|/gap
 896 *                    Note that both the residual and the gap are
 897 *                    proportional to the matrix, so ||T|| doesn't play
 898 *                    a role in the quotient
 899
 900 *
 901 *                    Convergence test for Rayleigh-Quotient iteration
 902 *                    (omitted when Bisection has been used)
 903 *
 904                      IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
 905      $                    RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
 906      $                    THEN
 907 *                       We need to check that the RQCORR update doesn't
 908 *                       move the eigenvalue away from the desired one and
 909 *                       towards a neighbor. -> protection with bisection
 910                         IF(INDEIG.LE.NEGCNT) THEN
 911 *                          The wanted eigenvalue lies to the left
 912                            SGNDEF = -ONE
 913                         ELSE
 914 *                          The wanted eigenvalue lies to the right
 915                            SGNDEF = ONE
 916                         ENDIF
 917 *                       We only use the RQCORR if it improves the
 918 *                       the iterate reasonably.
 919                         IF( ( RQCORR*SGNDEF.GE.ZERO )
 920      $                       .AND.( LAMBDA + RQCORR.LE. RIGHT)
 921      $                       .AND.( LAMBDA + RQCORR.GE. LEFT)
 922      $                       ) THEN
 923                            USEDRQ = .TRUE.
 924 *                          Store new midpoint of bisection interval in WORK
 925                            IF(SGNDEF.EQ.ONE) THEN
 926 *                             The current LAMBDA is on the left of the true
 927 *                             eigenvalue
 928                               LEFT = LAMBDA
 929 *                             We prefer to assume that the error estimate
 930 *                             is correct. We could make the interval not
 931 *                             as a bracket but to be modified if the RQCORR
 932 *                             chooses to. In this case, the RIGHT side should
 933 *                             be modified as follows:
 934 *                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
 935                            ELSE
 936 *                             The current LAMBDA is on the right of the true
 937 *                             eigenvalue
 938                               RIGHT = LAMBDA
 939 *                             See comment about assuming the error estimate is
 940 *                             correct above.
 941 *                              LEFT = MIN(LEFT, LAMBDA + RQCORR)
 942                            ENDIF
 943                            WORK( WINDEX ) =
 944      $                       HALF * (RIGHT + LEFT)
 945 *                          Take RQCORR since it has the correct sign and
 946 *                          improves the iterate reasonably
 947                            LAMBDA = LAMBDA + RQCORR
 948 *                          Update width of error interval
 949                            WERR( WINDEX ) =
 950      $                             HALF * (RIGHT-LEFT)
 951                         ELSE
 952                            NEEDBS = .TRUE.
 953                         ENDIF
 954                         IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
 955 *                             The eigenvalue is computed to bisection accuracy
 956 *                             compute eigenvector and stop
 957                            USEDBS = .TRUE.
 958                            GOTO 120
 959                         ELSEIF( ITER.LT.MAXITR ) THEN
 960                            GOTO 120
 961                         ELSEIF( ITER.EQ.MAXITR ) THEN
 962                            NEEDBS = .TRUE.
 963                            GOTO 120
 964                         ELSE
 965                            INFO = 5
 966                            RETURN
 967                         END IF
 968                      ELSE
 969                         STP2II = .FALSE.
 970         IF(USEDRQ .AND. USEDBS .AND.
 971      $                     BSTRES.LE.RESID) THEN
 972                            LAMBDA = BSTW
 973                            STP2II = .TRUE.
 974                         ENDIF
 975                         IF (STP2II) THEN
 976 *                          improve error angle by second step
 977                            CALL ZLAR1V( IN, 1, IN, LAMBDA,
 978      $                          D( IBEGIN ), L( IBEGIN ),
 979      $                          WORK(INDLD+IBEGIN-1),
 980      $                          WORK(INDLLD+IBEGIN-1),
 981      $                          PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
 982      $                          .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
 983      $                          IWORK( IINDR+WINDEX ),
 984      $                          ISUPPZ( 2*WINDEX-1 ),
 985      $                          NRMINV, RESID, RQCORR, WORK( INDWRK ) )
 986                         ENDIF
 987                         WORK( WINDEX ) = LAMBDA
 988                      END IF
 989 *
 990 *                    Compute FP-vector support w.r.t. whole matrix
 991 *
 992                      ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
 993                      ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
 994                      ZFROM = ISUPPZ( 2*WINDEX-1 )
 995                      ZTO = ISUPPZ( 2*WINDEX )
 996                      ISUPMN = ISUPMN + OLDIEN
 997                      ISUPMX = ISUPMX + OLDIEN
 998 *                    Ensure vector is ok if support in the RQI has changed
 999                      IF(ISUPMN.LT.ZFROM) THEN
1000                         DO 122 II = ISUPMN,ZFROM-1
1001                            Z( II, WINDEX ) = ZERO
1002  122                    CONTINUE
1003                      ENDIF
1004                      IF(ISUPMX.GT.ZTO) THEN
1005                         DO 123 II = ZTO+1,ISUPMX
1006                            Z( II, WINDEX ) = ZERO
1007  123                    CONTINUE
1008                      ENDIF
1009                      CALL ZDSCAL( ZTO-ZFROM+1, NRMINV,
1010      $                       Z( ZFROM, WINDEX ), 1 )
1011  125                 CONTINUE
1012 *                    Update W
1013                      W( WINDEX ) = LAMBDA+SIGMA
1014 *                    Recompute the gaps on the left and right
1015 *                    But only allow them to become larger and not
1016 *                    smaller (which can only happen through "bad"
1017 *                    cancellation and doesn't reflect the theory
1018 *                    where the initial gaps are underestimated due
1019 *                    to WERR being too crude.)
1020                      IF(.NOT.ESKIP) THEN
1021                         IF( K.GT.1) THEN
1022                            WGAP( WINDMN ) = MAX( WGAP(WINDMN),
1023      $                          W(WINDEX)-WERR(WINDEX)
1024      $                          - W(WINDMN)-WERR(WINDMN) )
1025                         ENDIF
1026                         IF( WINDEX.LT.WEND ) THEN
1027                            WGAP( WINDEX ) = MAX( SAVGAP,
1028      $                          W( WINDPL )-WERR( WINDPL )
1029      $                          - W( WINDEX )-WERR( WINDEX) )
1030                         ENDIF
1031                      ENDIF
1032                      IDONE = IDONE + 1
1033                   ENDIF
1034 *                 here ends the code for the current child
1035 *
1036  139              CONTINUE
1037 *                 Proceed to any remaining child nodes
1038                   NEWFST = J + 1
1039  140           CONTINUE
1040  150        CONTINUE
1041             NDEPTH = NDEPTH + 1
1042             GO TO 40
1043          END IF
1044          IBEGIN = IEND + 1
1045          WBEGIN = WEND + 1
1046  170  CONTINUE
1047 *
1048
1049       RETURN
1050 *
1051 *     End of ZLARRV
1052 *
1053       END