zlarrv.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 881 行 · 第 1/5 页

</span>                        NEWFTT = DOL - 1
                     ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
<span class="comment">*</span><span class="comment">                       Store representation at the right end of Z array
</span>                        NEWFTT = DOU
                     ELSE
                        NEWFTT = WBEGIN + NEWFST - 1
                     ENDIF
                  ENDIF

                  IF( NEWSIZ.GT.1) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Current child is not a singleton but a cluster.
</span><span class="comment">*</span><span class="comment">                    Compute and store new representation of child.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute left and right cluster gap.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    LGAP and RGAP are not computed from WORK because
</span><span class="comment">*</span><span class="comment">                    the eigenvalue approximations may stem from RRRs
</span><span class="comment">*</span><span class="comment">                    different shifts. However, W hold all eigenvalues
</span><span class="comment">*</span><span class="comment">                    of the unshifted matrix. Still, the entries in WGAP
</span><span class="comment">*</span><span class="comment">                    have to be computed from WORK since the entries
</span><span class="comment">*</span><span class="comment">                    in W might be of the same order so that gaps are not
</span><span class="comment">*</span><span class="comment">                    exhibited correctly for very close eigenvalues.
</span>                     IF( NEWFST.EQ.1 ) THEN
                        LGAP = MAX( ZERO,
     $                       W(WBEGIN)-WERR(WBEGIN) - VL )
                    ELSE
                        LGAP = WGAP( WBEGIN+NEWFST-2 )
                     ENDIF
                     RGAP = WGAP( WBEGIN+NEWLST-1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute left- and rightmost eigenvalue of child
</span><span class="comment">*</span><span class="comment">                    to high precision in order to shift as close
</span><span class="comment">*</span><span class="comment">                    as possible and obtain as large relative gaps
</span><span class="comment">*</span><span class="comment">                    as possible
</span><span class="comment">*</span><span class="comment">
</span>                     DO 55 K =1,2
                        IF(K.EQ.1) THEN
                           P = INDEXW( WBEGIN-1+NEWFST )
                        ELSE
                           P = INDEXW( WBEGIN-1+NEWLST )
                        ENDIF
                        OFFSET = INDEXW( WBEGIN ) - 1
                        CALL <a name="DLARRB.557"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>( IN, D(IBEGIN),
     $                       WORK( INDLLD+IBEGIN-1 ),P,P,
     $                       RQTOL, RQTOL, OFFSET,
     $                       WORK(WBEGIN),WGAP(WBEGIN),
     $                       WERR(WBEGIN),WORK( INDWRK ),
     $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
     $                       IN, IINFO )
 55                  CONTINUE
<span class="comment">*</span><span class="comment">
</span>                     IF((WBEGIN+NEWLST-1.LT.DOL).OR.
     $                  (WBEGIN+NEWFST-1.GT.DOU)) THEN
<span class="comment">*</span><span class="comment">                       if the cluster contains no desired eigenvalues
</span><span class="comment">*</span><span class="comment">                       skip the computation of that branch of the rep. tree
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       We could skip before the refinement of the extremal
</span><span class="comment">*</span><span class="comment">                       eigenvalues of the child, but then the representation
</span><span class="comment">*</span><span class="comment">                       tree could be different from the one when nothing is
</span><span class="comment">*</span><span class="comment">                       skipped. For this reason we skip at this place.
</span>                        IDONE = IDONE + NEWLST - NEWFST + 1
                        GOTO 139
                     ENDIF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute RRR of child cluster.
</span><span class="comment">*</span><span class="comment">                    Note that the new RRR is stored in Z
</span><span class="comment">*</span><span class="comment">
</span>C                    <a name="DLARRF.582"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a> needs LWORK = 2*N
                     CALL <a name="DLARRF.583"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>( IN, D( IBEGIN ), L( IBEGIN ),
     $                         WORK(INDLD+IBEGIN-1),
     $                         NEWFST, NEWLST, WORK(WBEGIN),
     $                         WGAP(WBEGIN), WERR(WBEGIN),
     $                         SPDIAM, LGAP, RGAP, PIVMIN, TAU,
     $                         WORK( INDIN1 ), WORK( INDIN2 ),
     $                         WORK( INDWRK ), IINFO )
<span class="comment">*</span><span class="comment">                    In the complex case, <a name="DLARRF.590"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a> cannot write
</span><span class="comment">*</span><span class="comment">                    the new RRR directly into Z and needs an intermediate
</span><span class="comment">*</span><span class="comment">                    workspace
</span>                     DO 56 K = 1, IN-1
                        Z( IBEGIN+K-1, NEWFTT ) =
     $                     DCMPLX( WORK( INDIN1+K-1 ), ZERO )
                        Z( IBEGIN+K-1, NEWFTT+1 ) =
     $                     DCMPLX( WORK( INDIN2+K-1 ), ZERO )
   56                CONTINUE
                     Z( IEND, NEWFTT ) =
     $                  DCMPLX( WORK( INDIN1+IN-1 ), ZERO )
                     IF( IINFO.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">                       a new RRR for the cluster was found by <a name="DLARRF.602"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>
</span><span class="comment">*</span><span class="comment">                       update shift and store it
</span>                        SSIGMA = SIGMA + TAU
                        Z( IEND, NEWFTT+1 ) = DCMPLX( SSIGMA, ZERO )
<span class="comment">*</span><span class="comment">                       WORK() are the midpoints and WERR() the semi-width
</span><span class="comment">*</span><span class="comment">                       Note that the entries in W are unchanged.
</span>                        DO 116 K = NEWFST, NEWLST
                           FUDGE =
     $                          THREE*EPS*ABS(WORK(WBEGIN+K-1))
                           WORK( WBEGIN + K - 1 ) =
     $                          WORK( WBEGIN + K - 1) - TAU
                           FUDGE = FUDGE +
     $                          FOUR*EPS*ABS(WORK(WBEGIN+K-1))
<span class="comment">*</span><span class="comment">                          Fudge errors
</span>                           WERR( WBEGIN + K - 1 ) =
     $                          WERR( WBEGIN + K - 1 ) + FUDGE
<span class="comment">*</span><span class="comment">                          Gaps are not fudged. Provided that WERR is small
</span><span class="comment">*</span><span class="comment">                          when eigenvalues are close, a zero gap indicates
</span><span class="comment">*</span><span class="comment">                          that a new representation is needed for resolving
</span><span class="comment">*</span><span class="comment">                          the cluster. A fudge could lead to a wrong decision
</span><span class="comment">*</span><span class="comment">                          of judging eigenvalues 'separated' which in
</span><span class="comment">*</span><span class="comment">                          reality are not. This could have a negative impact
</span><span class="comment">*</span><span class="comment">                          on the orthogonality of the computed eigenvectors.
</span> 116                    CONTINUE

                        NCLUS = NCLUS + 1
                        K = NEWCLS + 2*NCLUS
                        IWORK( K-1 ) = NEWFST
                        IWORK( K ) = NEWLST
                     ELSE
                        INFO = -2
                        RETURN
                     ENDIF
                  ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute eigenvector of singleton
</span><span class="comment">*</span><span class="comment">
</span>                     ITER = 0
<span class="comment">*</span><span class="comment">
</span>                     TOL = FOUR * LOG(DBLE(IN)) * EPS
<span class="comment">*</span><span class="comment">
</span>                     K = NEWFST
                     WINDEX = WBEGIN + K - 1
                     WINDMN = MAX(WINDEX - 1,1)
                     WINDPL = MIN(WINDEX + 1,M)
                     LAMBDA = WORK( WINDEX )
                     DONE = DONE + 1
<span class="comment">*</span><span class="comment">                    Check if eigenvector computation is to be skipped
</span>                     IF((WINDEX.LT.DOL).OR.
     $                  (WINDEX.GT.DOU)) THEN
                        ESKIP = .TRUE.
                        GOTO 125
                     ELSE
                        ESKIP = .FALSE.
                     ENDIF
                     LEFT = WORK( WINDEX ) - WERR( WINDEX )
                     RIGHT = WORK( WINDEX ) + WERR( WINDEX )
                     INDEIG = INDEXW( WINDEX )
<span class="comment">*</span><span class="comment">                    Note that since we compute the eigenpairs for a child,
</span><span class="comment">*</span><span class="comment">                    all eigenvalue approximations are w.r.t the same shift.
</span><span class="comment">*</span><span class="comment">                    In this case, the entries in WORK should be used for
</span><span class="comment">*</span><span class="comment">                    computing the gaps since they exhibit even very small
</span><span class="comment">*</span><span class="comment">                    differences in the eigenvalues, as opposed to the
</span><span class="comment">*</span><span class="comment">                    entries in W which might &quot;look&quot; the same.
</span>
                     IF( K .EQ. 1) THEN
<span class="comment">*</span><span class="comment">                       In the case RANGE='I' and with not much initial
</span><span class="comment">*</span><span class="comment">                       accuracy in LAMBDA and VL, the formula
</span><span class="comment">*</span><span class="comment">                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
</span><span class="comment">*</span><span class="comment">                       can lead to an overestimation of the left gap and
</span><span class="comment">*</span><span class="comment">                       thus to inadequately early RQI 'convergence'.
</span><span class="comment">*</span><span class="comment">                       Prevent this by forcing a small left gap.
</span>                        LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
                     ELSE
                        LGAP = WGAP(WINDMN)
                     ENDIF
                     IF( K .EQ. IM) THEN
<span class="comment">*</span><span class="comment">                       In the case RANGE='I' and with not much initial
</span><span class="comment">*</span><span class="comment">                       accuracy in LAMBDA and VU, the formula
</span><span class="comment">*</span><span class="comment">                       can lead to an overestimation of the right gap and
</span><span class="comment">*</span><span class="comment">                       thus to inadequately early RQI 'convergence'.
</span><span class="comment">*</span><span class="comment">                       Prevent this by forcing a small right gap.
</span>                        RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
                     ELSE
                        RGAP = WGAP(WINDEX)
                     ENDIF
                     GAP = MIN( LGAP, RGAP )
                     IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN