LAPACK 3.3.0

clarrv.f

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