LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlarrv.f
Go to the documentation of this file.
00001 *> \brief \b ZLARRV
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZLARRV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarrv.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarrv.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarrv.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZLARRV( N, VL, VU, D, L, PIVMIN,
00022 *                          ISPLIT, M, DOL, DOU, MINRGP,
00023 *                          RTOL1, RTOL2, W, WERR, WGAP,
00024 *                          IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
00025 *                          WORK, IWORK, INFO )
00026 * 
00027 *       .. Scalar Arguments ..
00028 *       INTEGER            DOL, DOU, INFO, LDZ, M, N
00029 *       DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
00033 *      $                   ISUPPZ( * ), IWORK( * )
00034 *       DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
00035 *      $                   WGAP( * ), WORK( * )
00036 *       COMPLEX*16        Z( LDZ, * )
00037 *       ..
00038 *  
00039 *
00040 *> \par Purpose:
00041 *  =============
00042 *>
00043 *> \verbatim
00044 *>
00045 *> ZLARRV computes the eigenvectors of the tridiagonal matrix
00046 *> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
00047 *> The input eigenvalues should have been computed by DLARRE.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] N
00054 *> \verbatim
00055 *>          N is INTEGER
00056 *>          The order of the matrix.  N >= 0.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] VL
00060 *> \verbatim
00061 *>          VL is DOUBLE PRECISION
00062 *> \endverbatim
00063 *>
00064 *> \param[in] VU
00065 *> \verbatim
00066 *>          VU is DOUBLE PRECISION
00067 *>          Lower and upper bounds of the interval that contains the desired
00068 *>          eigenvalues. VL < VU. Needed to compute gaps on the left or right
00069 *>          end of the extremal eigenvalues in the desired RANGE.
00070 *> \endverbatim
00071 *>
00072 *> \param[in,out] D
00073 *> \verbatim
00074 *>          D is DOUBLE PRECISION array, dimension (N)
00075 *>          On entry, the N diagonal elements of the diagonal matrix D.
00076 *>          On exit, D may be overwritten.
00077 *> \endverbatim
00078 *>
00079 *> \param[in,out] L
00080 *> \verbatim
00081 *>          L is DOUBLE PRECISION array, dimension (N)
00082 *>          On entry, the (N-1) subdiagonal elements of the unit
00083 *>          bidiagonal matrix L are in elements 1 to N-1 of L
00084 *>          (if the matrix is not splitted.) At the end of each block
00085 *>          is stored the corresponding shift as given by DLARRE.
00086 *>          On exit, L is overwritten.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] PIVMIN
00090 *> \verbatim
00091 *>          PIVMIN is DOUBLE PRECISION
00092 *>          The minimum pivot allowed in the Sturm sequence.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] ISPLIT
00096 *> \verbatim
00097 *>          ISPLIT is INTEGER array, dimension (N)
00098 *>          The splitting points, at which T breaks up into blocks.
00099 *>          The first block consists of rows/columns 1 to
00100 *>          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
00101 *>          through ISPLIT( 2 ), etc.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] M
00105 *> \verbatim
00106 *>          M is INTEGER
00107 *>          The total number of input eigenvalues.  0 <= M <= N.
00108 *> \endverbatim
00109 *>
00110 *> \param[in] DOL
00111 *> \verbatim
00112 *>          DOL is INTEGER
00113 *> \endverbatim
00114 *>
00115 *> \param[in] DOU
00116 *> \verbatim
00117 *>          DOU is INTEGER
00118 *>          If the user wants to compute only selected eigenvectors from all
00119 *>          the eigenvalues supplied, he can specify an index range DOL:DOU.
00120 *>          Or else the setting DOL=1, DOU=M should be applied.
00121 *>          Note that DOL and DOU refer to the order in which the eigenvalues
00122 *>          are stored in W.
00123 *>          If the user wants to compute only selected eigenpairs, then
00124 *>          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
00125 *>          computed eigenvectors. All other columns of Z are set to zero.
00126 *> \endverbatim
00127 *>
00128 *> \param[in] MINRGP
00129 *> \verbatim
00130 *>          MINRGP is DOUBLE PRECISION
00131 *> \endverbatim
00132 *>
00133 *> \param[in] RTOL1
00134 *> \verbatim
00135 *>          RTOL1 is DOUBLE PRECISION
00136 *> \endverbatim
00137 *>
00138 *> \param[in] RTOL2
00139 *> \verbatim
00140 *>          RTOL2 is DOUBLE PRECISION
00141 *>           Parameters for bisection.
00142 *>           An interval [LEFT,RIGHT] has converged if
00143 *>           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
00144 *> \endverbatim
00145 *>
00146 *> \param[in,out] W
00147 *> \verbatim
00148 *>          W is DOUBLE PRECISION array, dimension (N)
00149 *>          The first M elements of W contain the APPROXIMATE eigenvalues for
00150 *>          which eigenvectors are to be computed.  The eigenvalues
00151 *>          should be grouped by split-off block and ordered from
00152 *>          smallest to largest within the block ( The output array
00153 *>          W from DLARRE is expected here ). Furthermore, they are with
00154 *>          respect to the shift of the corresponding root representation
00155 *>          for their block. On exit, W holds the eigenvalues of the
00156 *>          UNshifted matrix.
00157 *> \endverbatim
00158 *>
00159 *> \param[in,out] WERR
00160 *> \verbatim
00161 *>          WERR is DOUBLE PRECISION array, dimension (N)
00162 *>          The first M elements contain the semiwidth of the uncertainty
00163 *>          interval of the corresponding eigenvalue in W
00164 *> \endverbatim
00165 *>
00166 *> \param[in,out] WGAP
00167 *> \verbatim
00168 *>          WGAP is DOUBLE PRECISION array, dimension (N)
00169 *>          The separation from the right neighbor eigenvalue in W.
00170 *> \endverbatim
00171 *>
00172 *> \param[in] IBLOCK
00173 *> \verbatim
00174 *>          IBLOCK is INTEGER array, dimension (N)
00175 *>          The indices of the blocks (submatrices) associated with the
00176 *>          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
00177 *>          W(i) belongs to the first block from the top, =2 if W(i)
00178 *>          belongs to the second block, etc.
00179 *> \endverbatim
00180 *>
00181 *> \param[in] INDEXW
00182 *> \verbatim
00183 *>          INDEXW is INTEGER array, dimension (N)
00184 *>          The indices of the eigenvalues within each block (submatrix);
00185 *>          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
00186 *>          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
00187 *> \endverbatim
00188 *>
00189 *> \param[in] GERS
00190 *> \verbatim
00191 *>          GERS is DOUBLE PRECISION array, dimension (2*N)
00192 *>          The N Gerschgorin intervals (the i-th Gerschgorin interval
00193 *>          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
00194 *>          be computed from the original UNshifted matrix.
00195 *> \endverbatim
00196 *>
00197 *> \param[out] Z
00198 *> \verbatim
00199 *>          Z is COMPLEX*16 array, dimension (LDZ, max(1,M) )
00200 *>          If INFO = 0, the first M columns of Z contain the
00201 *>          orthonormal eigenvectors of the matrix T
00202 *>          corresponding to the input eigenvalues, with the i-th
00203 *>          column of Z holding the eigenvector associated with W(i).
00204 *>          Note: the user must ensure that at least max(1,M) columns are
00205 *>          supplied in the array Z.
00206 *> \endverbatim
00207 *>
00208 *> \param[in] LDZ
00209 *> \verbatim
00210 *>          LDZ is INTEGER
00211 *>          The leading dimension of the array Z.  LDZ >= 1, and if
00212 *>          JOBZ = 'V', LDZ >= max(1,N).
00213 *> \endverbatim
00214 *>
00215 *> \param[out] ISUPPZ
00216 *> \verbatim
00217 *>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
00218 *>          The support of the eigenvectors in Z, i.e., the indices
00219 *>          indicating the nonzero elements in Z. The I-th eigenvector
00220 *>          is nonzero only in elements ISUPPZ( 2*I-1 ) through
00221 *>          ISUPPZ( 2*I ).
00222 *> \endverbatim
00223 *>
00224 *> \param[out] WORK
00225 *> \verbatim
00226 *>          WORK is DOUBLE PRECISION array, dimension (12*N)
00227 *> \endverbatim
00228 *>
00229 *> \param[out] IWORK
00230 *> \verbatim
00231 *>          IWORK is INTEGER array, dimension (7*N)
00232 *> \endverbatim
00233 *>
00234 *> \param[out] INFO
00235 *> \verbatim
00236 *>          INFO is INTEGER
00237 *>          = 0:  successful exit
00238 *>
00239 *>          > 0:  A problem occured in ZLARRV.
00240 *>          < 0:  One of the called subroutines signaled an internal problem.
00241 *>                Needs inspection of the corresponding parameter IINFO
00242 *>                for further information.
00243 *>
00244 *>          =-1:  Problem in DLARRB when refining a child's eigenvalues.
00245 *>          =-2:  Problem in DLARRF when computing the RRR of a child.
00246 *>                When a child is inside a tight cluster, it can be difficult
00247 *>                to find an RRR. A partial remedy from the user's point of
00248 *>                view is to make the parameter MINRGP smaller and recompile.
00249 *>                However, as the orthogonality of the computed vectors is
00250 *>                proportional to 1/MINRGP, the user should be aware that
00251 *>                he might be trading in precision when he decreases MINRGP.
00252 *>          =-3:  Problem in DLARRB when refining a single eigenvalue
00253 *>                after the Rayleigh correction was rejected.
00254 *>          = 5:  The Rayleigh Quotient Iteration failed to converge to
00255 *>                full accuracy in MAXITR steps.
00256 *> \endverbatim
00257 *
00258 *  Authors:
00259 *  ========
00260 *
00261 *> \author Univ. of Tennessee 
00262 *> \author Univ. of California Berkeley 
00263 *> \author Univ. of Colorado Denver 
00264 *> \author NAG Ltd. 
00265 *
00266 *> \date November 2011
00267 *
00268 *> \ingroup complex16OTHERauxiliary
00269 *
00270 *> \par Contributors:
00271 *  ==================
00272 *>
00273 *> Beresford Parlett, University of California, Berkeley, USA \n
00274 *> Jim Demmel, University of California, Berkeley, USA \n
00275 *> Inderjit Dhillon, University of Texas, Austin, USA \n
00276 *> Osni Marques, LBNL/NERSC, USA \n
00277 *> Christof Voemel, University of California, Berkeley, USA
00278 *
00279 *  =====================================================================
00280       SUBROUTINE ZLARRV( N, VL, VU, D, L, PIVMIN,
00281      $                   ISPLIT, M, DOL, DOU, MINRGP,
00282      $                   RTOL1, RTOL2, W, WERR, WGAP,
00283      $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
00284      $                   WORK, IWORK, INFO )
00285 *
00286 *  -- LAPACK auxiliary routine (version 3.4.0) --
00287 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00288 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00289 *     November 2011
00290 *
00291 *     .. Scalar Arguments ..
00292       INTEGER            DOL, DOU, INFO, LDZ, M, N
00293       DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
00294 *     ..
00295 *     .. Array Arguments ..
00296       INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
00297      $                   ISUPPZ( * ), IWORK( * )
00298       DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
00299      $                   WGAP( * ), WORK( * )
00300       COMPLEX*16        Z( LDZ, * )
00301 *     ..
00302 *
00303 *  =====================================================================
00304 *
00305 *     .. Parameters ..
00306       INTEGER            MAXITR
00307       PARAMETER          ( MAXITR = 10 )
00308       COMPLEX*16         CZERO
00309       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
00310       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, HALF
00311       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
00312      $                     TWO = 2.0D0, THREE = 3.0D0,
00313      $                     FOUR = 4.0D0, HALF = 0.5D0)
00314 *     ..
00315 *     .. Local Scalars ..
00316       LOGICAL            ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
00317       INTEGER            DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
00318      $                   IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
00319      $                   INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
00320      $                   ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
00321      $                   NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
00322      $                   NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
00323      $                   OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
00324      $                   WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
00325      $                   ZUSEDW
00326       INTEGER            INDIN1, INDIN2
00327       DOUBLE PRECISION   BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
00328      $                   LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
00329      $                   RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
00330      $                   SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
00331 *     ..
00332 *     .. External Functions ..
00333       DOUBLE PRECISION   DLAMCH
00334       EXTERNAL           DLAMCH
00335 *     ..
00336 *     .. External Subroutines ..
00337       EXTERNAL           DCOPY, DLARRB, DLARRF, ZDSCAL, ZLAR1V,
00338      $                   ZLASET
00339 *     ..
00340 *     .. Intrinsic Functions ..
00341       INTRINSIC ABS, DBLE, MAX, MIN
00342       INTRINSIC DCMPLX
00343 *     ..
00344 *     .. Executable Statements ..
00345 *     ..
00346 
00347 *     The first N entries of WORK are reserved for the eigenvalues
00348       INDLD = N+1
00349       INDLLD= 2*N+1
00350       INDIN1 = 3*N + 1
00351       INDIN2 = 4*N + 1
00352       INDWRK = 5*N + 1
00353       MINWSIZE = 12 * N
00354 
00355       DO 5 I= 1,MINWSIZE
00356          WORK( I ) = ZERO
00357  5    CONTINUE
00358 
00359 *     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
00360 *     factorization used to compute the FP vector
00361       IINDR = 0
00362 *     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
00363 *     layer and the one above.
00364       IINDC1 = N
00365       IINDC2 = 2*N
00366       IINDWK = 3*N + 1
00367 
00368       MINIWSIZE = 7 * N
00369       DO 10 I= 1,MINIWSIZE
00370          IWORK( I ) = 0
00371  10   CONTINUE
00372 
00373       ZUSEDL = 1
00374       IF(DOL.GT.1) THEN
00375 *        Set lower bound for use of Z
00376          ZUSEDL = DOL-1
00377       ENDIF
00378       ZUSEDU = M
00379       IF(DOU.LT.M) THEN
00380 *        Set lower bound for use of Z
00381          ZUSEDU = DOU+1
00382       ENDIF
00383 *     The width of the part of Z that is used
00384       ZUSEDW = ZUSEDU - ZUSEDL + 1
00385 
00386 
00387       CALL ZLASET( 'Full', N, ZUSEDW, CZERO, CZERO,
00388      $                    Z(1,ZUSEDL), LDZ )
00389 
00390       EPS = DLAMCH( 'Precision' )
00391       RQTOL = TWO * EPS
00392 *
00393 *     Set expert flags for standard code.
00394       TRYRQC = .TRUE.
00395 
00396       IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
00397       ELSE
00398 *        Only selected eigenpairs are computed. Since the other evalues
00399 *        are not refined by RQ iteration, bisection has to compute to full
00400 *        accuracy.
00401          RTOL1 = FOUR * EPS
00402          RTOL2 = FOUR * EPS
00403       ENDIF
00404 
00405 *     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
00406 *     desired eigenvalues. The support of the nonzero eigenvector
00407 *     entries is contained in the interval IBEGIN:IEND.
00408 *     Remark that if k eigenpairs are desired, then the eigenvectors
00409 *     are stored in k contiguous columns of Z.
00410 
00411 *     DONE is the number of eigenvectors already computed
00412       DONE = 0
00413       IBEGIN = 1
00414       WBEGIN = 1
00415       DO 170 JBLK = 1, IBLOCK( M )
00416          IEND = ISPLIT( JBLK )
00417          SIGMA = L( IEND )
00418 *        Find the eigenvectors of the submatrix indexed IBEGIN
00419 *        through IEND.
00420          WEND = WBEGIN - 1
00421  15      CONTINUE
00422          IF( WEND.LT.M ) THEN
00423             IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
00424                WEND = WEND + 1
00425                GO TO 15
00426             END IF
00427          END IF
00428          IF( WEND.LT.WBEGIN ) THEN
00429             IBEGIN = IEND + 1
00430             GO TO 170
00431          ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
00432             IBEGIN = IEND + 1
00433             WBEGIN = WEND + 1
00434             GO TO 170
00435          END IF
00436 
00437 *        Find local spectral diameter of the block
00438          GL = GERS( 2*IBEGIN-1 )
00439          GU = GERS( 2*IBEGIN )
00440          DO 20 I = IBEGIN+1 , IEND
00441             GL = MIN( GERS( 2*I-1 ), GL )
00442             GU = MAX( GERS( 2*I ), GU )
00443  20      CONTINUE
00444          SPDIAM = GU - GL
00445 
00446 *        OLDIEN is the last index of the previous block
00447          OLDIEN = IBEGIN - 1
00448 *        Calculate the size of the current block
00449          IN = IEND - IBEGIN + 1
00450 *        The number of eigenvalues in the current block
00451          IM = WEND - WBEGIN + 1
00452 
00453 *        This is for a 1x1 block
00454          IF( IBEGIN.EQ.IEND ) THEN
00455             DONE = DONE+1
00456             Z( IBEGIN, WBEGIN ) = DCMPLX( ONE, ZERO )
00457             ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
00458             ISUPPZ( 2*WBEGIN ) = IBEGIN
00459             W( WBEGIN ) = W( WBEGIN ) + SIGMA
00460             WORK( WBEGIN ) = W( WBEGIN )
00461             IBEGIN = IEND + 1
00462             WBEGIN = WBEGIN + 1
00463             GO TO 170
00464          END IF
00465 
00466 *        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
00467 *        Note that these can be approximations, in this case, the corresp.
00468 *        entries of WERR give the size of the uncertainty interval.
00469 *        The eigenvalue approximations will be refined when necessary as
00470 *        high relative accuracy is required for the computation of the
00471 *        corresponding eigenvectors.
00472          CALL DCOPY( IM, W( WBEGIN ), 1,
00473      $                   WORK( WBEGIN ), 1 )
00474 
00475 *        We store in W the eigenvalue approximations w.r.t. the original
00476 *        matrix T.
00477          DO 30 I=1,IM
00478             W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
00479  30      CONTINUE
00480 
00481 
00482 *        NDEPTH is the current depth of the representation tree
00483          NDEPTH = 0
00484 *        PARITY is either 1 or 0
00485          PARITY = 1
00486 *        NCLUS is the number of clusters for the next level of the
00487 *        representation tree, we start with NCLUS = 1 for the root
00488          NCLUS = 1
00489          IWORK( IINDC1+1 ) = 1
00490          IWORK( IINDC1+2 ) = IM
00491 
00492 *        IDONE is the number of eigenvectors already computed in the current
00493 *        block
00494          IDONE = 0
00495 *        loop while( IDONE.LT.IM )
00496 *        generate the representation tree for the current block and
00497 *        compute the eigenvectors
00498    40    CONTINUE
00499          IF( IDONE.LT.IM ) THEN
00500 *           This is a crude protection against infinitely deep trees
00501             IF( NDEPTH.GT.M ) THEN
00502                INFO = -2
00503                RETURN
00504             ENDIF
00505 *           breadth first processing of the current level of the representation
00506 *           tree: OLDNCL = number of clusters on current level
00507             OLDNCL = NCLUS
00508 *           reset NCLUS to count the number of child clusters
00509             NCLUS = 0
00510 *
00511             PARITY = 1 - PARITY
00512             IF( PARITY.EQ.0 ) THEN
00513                OLDCLS = IINDC1
00514                NEWCLS = IINDC2
00515             ELSE
00516                OLDCLS = IINDC2
00517                NEWCLS = IINDC1
00518             END IF
00519 *           Process the clusters on the current level
00520             DO 150 I = 1, OLDNCL
00521                J = OLDCLS + 2*I
00522 *              OLDFST, OLDLST = first, last index of current cluster.
00523 *                               cluster indices start with 1 and are relative
00524 *                               to WBEGIN when accessing W, WGAP, WERR, Z
00525                OLDFST = IWORK( J-1 )
00526                OLDLST = IWORK( J )
00527                IF( NDEPTH.GT.0 ) THEN
00528 *                 Retrieve relatively robust representation (RRR) of cluster
00529 *                 that has been computed at the previous level
00530 *                 The RRR is stored in Z and overwritten once the eigenvectors
00531 *                 have been computed or when the cluster is refined
00532 
00533                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
00534 *                    Get representation from location of the leftmost evalue
00535 *                    of the cluster
00536                      J = WBEGIN + OLDFST - 1
00537                   ELSE
00538                      IF(WBEGIN+OLDFST-1.LT.DOL) THEN
00539 *                       Get representation from the left end of Z array
00540                         J = DOL - 1
00541                      ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
00542 *                       Get representation from the right end of Z array
00543                         J = DOU
00544                      ELSE
00545                         J = WBEGIN + OLDFST - 1
00546                      ENDIF
00547                   ENDIF
00548                   DO 45 K = 1, IN - 1
00549                      D( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1,
00550      $                                 J ) )
00551                      L( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1,
00552      $                                 J+1 ) )
00553    45             CONTINUE
00554                   D( IEND ) = DBLE( Z( IEND, J ) )
00555                   SIGMA = DBLE( Z( IEND, J+1 ) )
00556 
00557 *                 Set the corresponding entries in Z to zero
00558                   CALL ZLASET( 'Full', IN, 2, CZERO, CZERO,
00559      $                         Z( IBEGIN, J), LDZ )
00560                END IF
00561 
00562 *              Compute DL and DLL of current RRR
00563                DO 50 J = IBEGIN, IEND-1
00564                   TMP = D( J )*L( J )
00565                   WORK( INDLD-1+J ) = TMP
00566                   WORK( INDLLD-1+J ) = TMP*L( J )
00567    50          CONTINUE
00568 
00569                IF( NDEPTH.GT.0 ) THEN
00570 *                 P and Q are index of the first and last eigenvalue to compute
00571 *                 within the current block
00572                   P = INDEXW( WBEGIN-1+OLDFST )
00573                   Q = INDEXW( WBEGIN-1+OLDLST )
00574 *                 Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
00575 *                 through the Q-OFFSET elements of these arrays are to be used.
00576 *                  OFFSET = P-OLDFST
00577                   OFFSET = INDEXW( WBEGIN ) - 1
00578 *                 perform limited bisection (if necessary) to get approximate
00579 *                 eigenvalues to the precision needed.
00580                   CALL DLARRB( IN, D( IBEGIN ),
00581      $                         WORK(INDLLD+IBEGIN-1),
00582      $                         P, Q, RTOL1, RTOL2, OFFSET,
00583      $                         WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
00584      $                         WORK( INDWRK ), IWORK( IINDWK ),
00585      $                         PIVMIN, SPDIAM, IN, IINFO )
00586                   IF( IINFO.NE.0 ) THEN
00587                      INFO = -1
00588                      RETURN
00589                   ENDIF
00590 *                 We also recompute the extremal gaps. W holds all eigenvalues
00591 *                 of the unshifted matrix and must be used for computation
00592 *                 of WGAP, the entries of WORK might stem from RRRs with
00593 *                 different shifts. The gaps from WBEGIN-1+OLDFST to
00594 *                 WBEGIN-1+OLDLST are correctly computed in DLARRB.
00595 *                 However, we only allow the gaps to become greater since
00596 *                 this is what should happen when we decrease WERR
00597                   IF( OLDFST.GT.1) THEN
00598                      WGAP( WBEGIN+OLDFST-2 ) =
00599      $             MAX(WGAP(WBEGIN+OLDFST-2),
00600      $                 W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
00601      $                 - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
00602                   ENDIF
00603                   IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
00604                      WGAP( WBEGIN+OLDLST-1 ) =
00605      $               MAX(WGAP(WBEGIN+OLDLST-1),
00606      $                   W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
00607      $                   - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
00608                   ENDIF
00609 *                 Each time the eigenvalues in WORK get refined, we store
00610 *                 the newly found approximation with all shifts applied in W
00611                   DO 53 J=OLDFST,OLDLST
00612                      W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
00613  53               CONTINUE
00614                END IF
00615 
00616 *              Process the current node.
00617                NEWFST = OLDFST
00618                DO 140 J = OLDFST, OLDLST
00619                   IF( J.EQ.OLDLST ) THEN
00620 *                    we are at the right end of the cluster, this is also the
00621 *                    boundary of the child cluster
00622                      NEWLST = J
00623                   ELSE IF ( WGAP( WBEGIN + J -1).GE.
00624      $                    MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
00625 *                    the right relative gap is big enough, the child cluster
00626 *                    (NEWFST,..,NEWLST) is well separated from the following
00627                      NEWLST = J
00628                    ELSE
00629 *                    inside a child cluster, the relative gap is not
00630 *                    big enough.
00631                      GOTO 140
00632                   END IF
00633 
00634 *                 Compute size of child cluster found
00635                   NEWSIZ = NEWLST - NEWFST + 1
00636 
00637 *                 NEWFTT is the place in Z where the new RRR or the computed
00638 *                 eigenvector is to be stored
00639                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
00640 *                    Store representation at location of the leftmost evalue
00641 *                    of the cluster
00642                      NEWFTT = WBEGIN + NEWFST - 1
00643                   ELSE
00644                      IF(WBEGIN+NEWFST-1.LT.DOL) THEN
00645 *                       Store representation at the left end of Z array
00646                         NEWFTT = DOL - 1
00647                      ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
00648 *                       Store representation at the right end of Z array
00649                         NEWFTT = DOU
00650                      ELSE
00651                         NEWFTT = WBEGIN + NEWFST - 1
00652                      ENDIF
00653                   ENDIF
00654 
00655                   IF( NEWSIZ.GT.1) THEN
00656 *
00657 *                    Current child is not a singleton but a cluster.
00658 *                    Compute and store new representation of child.
00659 *
00660 *
00661 *                    Compute left and right cluster gap.
00662 *
00663 *                    LGAP and RGAP are not computed from WORK because
00664 *                    the eigenvalue approximations may stem from RRRs
00665 *                    different shifts. However, W hold all eigenvalues
00666 *                    of the unshifted matrix. Still, the entries in WGAP
00667 *                    have to be computed from WORK since the entries
00668 *                    in W might be of the same order so that gaps are not
00669 *                    exhibited correctly for very close eigenvalues.
00670                      IF( NEWFST.EQ.1 ) THEN
00671                         LGAP = MAX( ZERO,
00672      $                       W(WBEGIN)-WERR(WBEGIN) - VL )
00673                     ELSE
00674                         LGAP = WGAP( WBEGIN+NEWFST-2 )
00675                      ENDIF
00676                      RGAP = WGAP( WBEGIN+NEWLST-1 )
00677 *
00678 *                    Compute left- and rightmost eigenvalue of child
00679 *                    to high precision in order to shift as close
00680 *                    as possible and obtain as large relative gaps
00681 *                    as possible
00682 *
00683                      DO 55 K =1,2
00684                         IF(K.EQ.1) THEN
00685                            P = INDEXW( WBEGIN-1+NEWFST )
00686                         ELSE
00687                            P = INDEXW( WBEGIN-1+NEWLST )
00688                         ENDIF
00689                         OFFSET = INDEXW( WBEGIN ) - 1
00690                         CALL DLARRB( IN, D(IBEGIN),
00691      $                       WORK( INDLLD+IBEGIN-1 ),P,P,
00692      $                       RQTOL, RQTOL, OFFSET,
00693      $                       WORK(WBEGIN),WGAP(WBEGIN),
00694      $                       WERR(WBEGIN),WORK( INDWRK ),
00695      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
00696      $                       IN, IINFO )
00697  55                  CONTINUE
00698 *
00699                      IF((WBEGIN+NEWLST-1.LT.DOL).OR.
00700      $                  (WBEGIN+NEWFST-1.GT.DOU)) THEN
00701 *                       if the cluster contains no desired eigenvalues
00702 *                       skip the computation of that branch of the rep. tree
00703 *
00704 *                       We could skip before the refinement of the extremal
00705 *                       eigenvalues of the child, but then the representation
00706 *                       tree could be different from the one when nothing is
00707 *                       skipped. For this reason we skip at this place.
00708                         IDONE = IDONE + NEWLST - NEWFST + 1
00709                         GOTO 139
00710                      ENDIF
00711 *
00712 *                    Compute RRR of child cluster.
00713 *                    Note that the new RRR is stored in Z
00714 *
00715 *                    DLARRF needs LWORK = 2*N
00716                      CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ),
00717      $                         WORK(INDLD+IBEGIN-1),
00718      $                         NEWFST, NEWLST, WORK(WBEGIN),
00719      $                         WGAP(WBEGIN), WERR(WBEGIN),
00720      $                         SPDIAM, LGAP, RGAP, PIVMIN, TAU,
00721      $                         WORK( INDIN1 ), WORK( INDIN2 ),
00722      $                         WORK( INDWRK ), IINFO )
00723 *                    In the complex case, DLARRF cannot write
00724 *                    the new RRR directly into Z and needs an intermediate
00725 *                    workspace
00726                      DO 56 K = 1, IN-1
00727                         Z( IBEGIN+K-1, NEWFTT ) =
00728      $                     DCMPLX( WORK( INDIN1+K-1 ), ZERO )
00729                         Z( IBEGIN+K-1, NEWFTT+1 ) =
00730      $                     DCMPLX( WORK( INDIN2+K-1 ), ZERO )
00731    56                CONTINUE
00732                      Z( IEND, NEWFTT ) =
00733      $                  DCMPLX( WORK( INDIN1+IN-1 ), ZERO )
00734                      IF( IINFO.EQ.0 ) THEN
00735 *                       a new RRR for the cluster was found by DLARRF
00736 *                       update shift and store it
00737                         SSIGMA = SIGMA + TAU
00738                         Z( IEND, NEWFTT+1 ) = DCMPLX( SSIGMA, ZERO )
00739 *                       WORK() are the midpoints and WERR() the semi-width
00740 *                       Note that the entries in W are unchanged.
00741                         DO 116 K = NEWFST, NEWLST
00742                            FUDGE =
00743      $                          THREE*EPS*ABS(WORK(WBEGIN+K-1))
00744                            WORK( WBEGIN + K - 1 ) =
00745      $                          WORK( WBEGIN + K - 1) - TAU
00746                            FUDGE = FUDGE +
00747      $                          FOUR*EPS*ABS(WORK(WBEGIN+K-1))
00748 *                          Fudge errors
00749                            WERR( WBEGIN + K - 1 ) =
00750      $                          WERR( WBEGIN + K - 1 ) + FUDGE
00751 *                          Gaps are not fudged. Provided that WERR is small
00752 *                          when eigenvalues are close, a zero gap indicates
00753 *                          that a new representation is needed for resolving
00754 *                          the cluster. A fudge could lead to a wrong decision
00755 *                          of judging eigenvalues 'separated' which in
00756 *                          reality are not. This could have a negative impact
00757 *                          on the orthogonality of the computed eigenvectors.
00758  116                    CONTINUE
00759 
00760                         NCLUS = NCLUS + 1
00761                         K = NEWCLS + 2*NCLUS
00762                         IWORK( K-1 ) = NEWFST
00763                         IWORK( K ) = NEWLST
00764                      ELSE
00765                         INFO = -2
00766                         RETURN
00767                      ENDIF
00768                   ELSE
00769 *
00770 *                    Compute eigenvector of singleton
00771 *
00772                      ITER = 0
00773 *
00774                      TOL = FOUR * LOG(DBLE(IN)) * EPS
00775 *
00776                      K = NEWFST
00777                      WINDEX = WBEGIN + K - 1
00778                      WINDMN = MAX(WINDEX - 1,1)
00779                      WINDPL = MIN(WINDEX + 1,M)
00780                      LAMBDA = WORK( WINDEX )
00781                      DONE = DONE + 1
00782 *                    Check if eigenvector computation is to be skipped
00783                      IF((WINDEX.LT.DOL).OR.
00784      $                  (WINDEX.GT.DOU)) THEN
00785                         ESKIP = .TRUE.
00786                         GOTO 125
00787                      ELSE
00788                         ESKIP = .FALSE.
00789                      ENDIF
00790                      LEFT = WORK( WINDEX ) - WERR( WINDEX )
00791                      RIGHT = WORK( WINDEX ) + WERR( WINDEX )
00792                      INDEIG = INDEXW( WINDEX )
00793 *                    Note that since we compute the eigenpairs for a child,
00794 *                    all eigenvalue approximations are w.r.t the same shift.
00795 *                    In this case, the entries in WORK should be used for
00796 *                    computing the gaps since they exhibit even very small
00797 *                    differences in the eigenvalues, as opposed to the
00798 *                    entries in W which might "look" the same.
00799 
00800                      IF( K .EQ. 1) THEN
00801 *                       In the case RANGE='I' and with not much initial
00802 *                       accuracy in LAMBDA and VL, the formula
00803 *                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
00804 *                       can lead to an overestimation of the left gap and
00805 *                       thus to inadequately early RQI 'convergence'.
00806 *                       Prevent this by forcing a small left gap.
00807                         LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
00808                      ELSE
00809                         LGAP = WGAP(WINDMN)
00810                      ENDIF
00811                      IF( K .EQ. IM) THEN
00812 *                       In the case RANGE='I' and with not much initial
00813 *                       accuracy in LAMBDA and VU, the formula
00814 *                       can lead to an overestimation of the right gap and
00815 *                       thus to inadequately early RQI 'convergence'.
00816 *                       Prevent this by forcing a small right gap.
00817                         RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
00818                      ELSE
00819                         RGAP = WGAP(WINDEX)
00820                      ENDIF
00821                      GAP = MIN( LGAP, RGAP )
00822                      IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN
00823 *                       The eigenvector support can become wrong
00824 *                       because significant entries could be cut off due to a
00825 *                       large GAPTOL parameter in LAR1V. Prevent this.
00826                         GAPTOL = ZERO
00827                      ELSE
00828                         GAPTOL = GAP * EPS
00829                      ENDIF
00830                      ISUPMN = IN
00831                      ISUPMX = 1
00832 *                    Update WGAP so that it holds the minimum gap
00833 *                    to the left or the right. This is crucial in the
00834 *                    case where bisection is used to ensure that the
00835 *                    eigenvalue is refined up to the required precision.
00836 *                    The correct value is restored afterwards.
00837                      SAVGAP = WGAP(WINDEX)
00838                      WGAP(WINDEX) = GAP
00839 *                    We want to use the Rayleigh Quotient Correction
00840 *                    as often as possible since it converges quadratically
00841 *                    when we are close enough to the desired eigenvalue.
00842 *                    However, the Rayleigh Quotient can have the wrong sign
00843 *                    and lead us away from the desired eigenvalue. In this
00844 *                    case, the best we can do is to use bisection.
00845                      USEDBS = .FALSE.
00846                      USEDRQ = .FALSE.
00847 *                    Bisection is initially turned off unless it is forced
00848                      NEEDBS =  .NOT.TRYRQC
00849  120                 CONTINUE
00850 *                    Check if bisection should be used to refine eigenvalue
00851                      IF(NEEDBS) THEN
00852 *                       Take the bisection as new iterate
00853                         USEDBS = .TRUE.
00854                         ITMP1 = IWORK( IINDR+WINDEX )
00855                         OFFSET = INDEXW( WBEGIN ) - 1
00856                         CALL DLARRB( IN, D(IBEGIN),
00857      $                       WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
00858      $                       ZERO, TWO*EPS, OFFSET,
00859      $                       WORK(WBEGIN),WGAP(WBEGIN),
00860      $                       WERR(WBEGIN),WORK( INDWRK ),
00861      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
00862      $                       ITMP1, IINFO )
00863                         IF( IINFO.NE.0 ) THEN
00864                            INFO = -3
00865                            RETURN
00866                         ENDIF
00867                         LAMBDA = WORK( WINDEX )
00868 *                       Reset twist index from inaccurate LAMBDA to
00869 *                       force computation of true MINGMA
00870                         IWORK( IINDR+WINDEX ) = 0
00871                      ENDIF
00872 *                    Given LAMBDA, compute the eigenvector.
00873                      CALL ZLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ),
00874      $                    L( IBEGIN ), WORK(INDLD+IBEGIN-1),
00875      $                    WORK(INDLLD+IBEGIN-1),
00876      $                    PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
00877      $                    .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
00878      $                    IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
00879      $                    NRMINV, RESID, RQCORR, WORK( INDWRK ) )
00880                      IF(ITER .EQ. 0) THEN
00881                         BSTRES = RESID
00882                         BSTW = LAMBDA
00883                      ELSEIF(RESID.LT.BSTRES) THEN
00884                         BSTRES = RESID
00885                         BSTW = LAMBDA
00886                      ENDIF
00887                      ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
00888                      ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
00889                      ITER = ITER + 1
00890 
00891 *                    sin alpha <= |resid|/gap
00892 *                    Note that both the residual and the gap are
00893 *                    proportional to the matrix, so ||T|| doesn't play
00894 *                    a role in the quotient
00895 
00896 *
00897 *                    Convergence test for Rayleigh-Quotient iteration
00898 *                    (omitted when Bisection has been used)
00899 *
00900                      IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
00901      $                    RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
00902      $                    THEN
00903 *                       We need to check that the RQCORR update doesn't
00904 *                       move the eigenvalue away from the desired one and
00905 *                       towards a neighbor. -> protection with bisection
00906                         IF(INDEIG.LE.NEGCNT) THEN
00907 *                          The wanted eigenvalue lies to the left
00908                            SGNDEF = -ONE
00909                         ELSE
00910 *                          The wanted eigenvalue lies to the right
00911                            SGNDEF = ONE
00912                         ENDIF
00913 *                       We only use the RQCORR if it improves the
00914 *                       the iterate reasonably.
00915                         IF( ( RQCORR*SGNDEF.GE.ZERO )
00916      $                       .AND.( LAMBDA + RQCORR.LE. RIGHT)
00917      $                       .AND.( LAMBDA + RQCORR.GE. LEFT)
00918      $                       ) THEN
00919                            USEDRQ = .TRUE.
00920 *                          Store new midpoint of bisection interval in WORK
00921                            IF(SGNDEF.EQ.ONE) THEN
00922 *                             The current LAMBDA is on the left of the true
00923 *                             eigenvalue
00924                               LEFT = LAMBDA
00925 *                             We prefer to assume that the error estimate
00926 *                             is correct. We could make the interval not
00927 *                             as a bracket but to be modified if the RQCORR
00928 *                             chooses to. In this case, the RIGHT side should
00929 *                             be modified as follows:
00930 *                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
00931                            ELSE
00932 *                             The current LAMBDA is on the right of the true
00933 *                             eigenvalue
00934                               RIGHT = LAMBDA
00935 *                             See comment about assuming the error estimate is
00936 *                             correct above.
00937 *                              LEFT = MIN(LEFT, LAMBDA + RQCORR)
00938                            ENDIF
00939                            WORK( WINDEX ) =
00940      $                       HALF * (RIGHT + LEFT)
00941 *                          Take RQCORR since it has the correct sign and
00942 *                          improves the iterate reasonably
00943                            LAMBDA = LAMBDA + RQCORR
00944 *                          Update width of error interval
00945                            WERR( WINDEX ) =
00946      $                             HALF * (RIGHT-LEFT)
00947                         ELSE
00948                            NEEDBS = .TRUE.
00949                         ENDIF
00950                         IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
00951 *                             The eigenvalue is computed to bisection accuracy
00952 *                             compute eigenvector and stop
00953                            USEDBS = .TRUE.
00954                            GOTO 120
00955                         ELSEIF( ITER.LT.MAXITR ) THEN
00956                            GOTO 120
00957                         ELSEIF( ITER.EQ.MAXITR ) THEN
00958                            NEEDBS = .TRUE.
00959                            GOTO 120
00960                         ELSE
00961                            INFO = 5
00962                            RETURN
00963                         END IF
00964                      ELSE
00965                         STP2II = .FALSE.
00966         IF(USEDRQ .AND. USEDBS .AND.
00967      $                     BSTRES.LE.RESID) THEN
00968                            LAMBDA = BSTW
00969                            STP2II = .TRUE.
00970                         ENDIF
00971                         IF (STP2II) THEN
00972 *                          improve error angle by second step
00973                            CALL ZLAR1V( IN, 1, IN, LAMBDA,
00974      $                          D( IBEGIN ), L( IBEGIN ),
00975      $                          WORK(INDLD+IBEGIN-1),
00976      $                          WORK(INDLLD+IBEGIN-1),
00977      $                          PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
00978      $                          .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
00979      $                          IWORK( IINDR+WINDEX ),
00980      $                          ISUPPZ( 2*WINDEX-1 ),
00981      $                          NRMINV, RESID, RQCORR, WORK( INDWRK ) )
00982                         ENDIF
00983                         WORK( WINDEX ) = LAMBDA
00984                      END IF
00985 *
00986 *                    Compute FP-vector support w.r.t. whole matrix
00987 *
00988                      ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
00989                      ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
00990                      ZFROM = ISUPPZ( 2*WINDEX-1 )
00991                      ZTO = ISUPPZ( 2*WINDEX )
00992                      ISUPMN = ISUPMN + OLDIEN
00993                      ISUPMX = ISUPMX + OLDIEN
00994 *                    Ensure vector is ok if support in the RQI has changed
00995                      IF(ISUPMN.LT.ZFROM) THEN
00996                         DO 122 II = ISUPMN,ZFROM-1
00997                            Z( II, WINDEX ) = ZERO
00998  122                    CONTINUE
00999                      ENDIF
01000                      IF(ISUPMX.GT.ZTO) THEN
01001                         DO 123 II = ZTO+1,ISUPMX
01002                            Z( II, WINDEX ) = ZERO
01003  123                    CONTINUE
01004                      ENDIF
01005                      CALL ZDSCAL( ZTO-ZFROM+1, NRMINV,
01006      $                       Z( ZFROM, WINDEX ), 1 )
01007  125                 CONTINUE
01008 *                    Update W
01009                      W( WINDEX ) = LAMBDA+SIGMA
01010 *                    Recompute the gaps on the left and right
01011 *                    But only allow them to become larger and not
01012 *                    smaller (which can only happen through "bad"
01013 *                    cancellation and doesn't reflect the theory
01014 *                    where the initial gaps are underestimated due
01015 *                    to WERR being too crude.)
01016                      IF(.NOT.ESKIP) THEN
01017                         IF( K.GT.1) THEN
01018                            WGAP( WINDMN ) = MAX( WGAP(WINDMN),
01019      $                          W(WINDEX)-WERR(WINDEX)
01020      $                          - W(WINDMN)-WERR(WINDMN) )
01021                         ENDIF
01022                         IF( WINDEX.LT.WEND ) THEN
01023                            WGAP( WINDEX ) = MAX( SAVGAP,
01024      $                          W( WINDPL )-WERR( WINDPL )
01025      $                          - W( WINDEX )-WERR( WINDEX) )
01026                         ENDIF
01027                      ENDIF
01028                      IDONE = IDONE + 1
01029                   ENDIF
01030 *                 here ends the code for the current child
01031 *
01032  139              CONTINUE
01033 *                 Proceed to any remaining child nodes
01034                   NEWFST = J + 1
01035  140           CONTINUE
01036  150        CONTINUE
01037             NDEPTH = NDEPTH + 1
01038             GO TO 40
01039          END IF
01040          IBEGIN = IEND + 1
01041          WBEGIN = WEND + 1
01042  170  CONTINUE
01043 *
01044 
01045       RETURN
01046 *
01047 *     End of ZLARRV
01048 *
01049       END
 All Files Functions