zlarrv.f.html

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

<span class="comment">*</span><span class="comment">                       The eigenvector support can become wrong
</span><span class="comment">*</span><span class="comment">                       because significant entries could be cut off due to a
</span><span class="comment">*</span><span class="comment">                       large GAPTOL parameter in LAR1V. Prevent this.
</span>                        GAPTOL = ZERO
                     ELSE
                        GAPTOL = GAP * EPS
                     ENDIF
                     ISUPMN = IN
                     ISUPMX = 1
<span class="comment">*</span><span class="comment">                    Update WGAP so that it holds the minimum gap
</span><span class="comment">*</span><span class="comment">                    to the left or the right. This is crucial in the
</span><span class="comment">*</span><span class="comment">                    case where bisection is used to ensure that the
</span><span class="comment">*</span><span class="comment">                    eigenvalue is refined up to the required precision.
</span><span class="comment">*</span><span class="comment">                    The correct value is restored afterwards.
</span>                     SAVGAP = WGAP(WINDEX)
                     WGAP(WINDEX) = GAP
<span class="comment">*</span><span class="comment">                    We want to use the Rayleigh Quotient Correction
</span><span class="comment">*</span><span class="comment">                    as often as possible since it converges quadratically
</span><span class="comment">*</span><span class="comment">                    when we are close enough to the desired eigenvalue.
</span><span class="comment">*</span><span class="comment">                    However, the Rayleigh Quotient can have the wrong sign
</span><span class="comment">*</span><span class="comment">                    and lead us away from the desired eigenvalue. In this
</span><span class="comment">*</span><span class="comment">                    case, the best we can do is to use bisection.
</span>                     USEDBS = .FALSE.
                     USEDRQ = .FALSE.
<span class="comment">*</span><span class="comment">                    Bisection is initially turned off unless it is forced
</span>                     NEEDBS =  .NOT.TRYRQC
 120                 CONTINUE
<span class="comment">*</span><span class="comment">                    Check if bisection should be used to refine eigenvalue
</span>                     IF(NEEDBS) THEN
<span class="comment">*</span><span class="comment">                       Take the bisection as new iterate
</span>                        USEDBS = .TRUE.
                        ITMP1 = IWORK( IINDR+WINDEX )
                        OFFSET = INDEXW( WBEGIN ) - 1
                        CALL <a name="DLARRB.723"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>( IN, D(IBEGIN),
     $                       WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
     $                       ZERO, TWO*EPS, OFFSET,
     $                       WORK(WBEGIN),WGAP(WBEGIN),
     $                       WERR(WBEGIN),WORK( INDWRK ),
     $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
     $                       ITMP1, IINFO )
                        IF( IINFO.NE.0 ) THEN
                           INFO = -3
                           RETURN
                        ENDIF
                        LAMBDA = WORK( WINDEX )
<span class="comment">*</span><span class="comment">                       Reset twist index from inaccurate LAMBDA to
</span><span class="comment">*</span><span class="comment">                       force computation of true MINGMA
</span>                        IWORK( IINDR+WINDEX ) = 0
                     ENDIF
<span class="comment">*</span><span class="comment">                    Given LAMBDA, compute the eigenvector.
</span>                     CALL <a name="ZLAR1V.740"></a><a href="zlar1v.f.html#ZLAR1V.1">ZLAR1V</a>( IN, 1, IN, LAMBDA, D( IBEGIN ),
     $                    L( IBEGIN ), WORK(INDLD+IBEGIN-1),
     $                    WORK(INDLLD+IBEGIN-1),
     $                    PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
     $                    .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
     $                    IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
     $                    NRMINV, RESID, RQCORR, WORK( INDWRK ) )
                     IF(ITER .EQ. 0) THEN
                        BSTRES = RESID
                        BSTW = LAMBDA
                     ELSEIF(RESID.LT.BSTRES) THEN
                        BSTRES = RESID
                        BSTW = LAMBDA
                     ENDIF
                     ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
                     ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
                     ITER = ITER + 1

<span class="comment">*</span><span class="comment">                    sin alpha &lt;= |resid|/gap
</span><span class="comment">*</span><span class="comment">                    Note that both the residual and the gap are
</span><span class="comment">*</span><span class="comment">                    proportional to the matrix, so ||T|| doesn't play
</span><span class="comment">*</span><span class="comment">                    a role in the quotient
</span>
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Convergence test for Rayleigh-Quotient iteration
</span><span class="comment">*</span><span class="comment">                    (omitted when Bisection has been used)
</span><span class="comment">*</span><span class="comment">
</span>                     IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
     $                    RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
     $                    THEN
<span class="comment">*</span><span class="comment">                       We need to check that the RQCORR update doesn't
</span><span class="comment">*</span><span class="comment">                       move the eigenvalue away from the desired one and
</span><span class="comment">*</span><span class="comment">                       towards a neighbor. -&gt; protection with bisection
</span>                        IF(INDEIG.LE.NEGCNT) THEN
<span class="comment">*</span><span class="comment">                          The wanted eigenvalue lies to the left
</span>                           SGNDEF = -ONE
                        ELSE
<span class="comment">*</span><span class="comment">                          The wanted eigenvalue lies to the right
</span>                           SGNDEF = ONE
                        ENDIF
<span class="comment">*</span><span class="comment">                       We only use the RQCORR if it improves the
</span><span class="comment">*</span><span class="comment">                       the iterate reasonably.
</span>                        IF( ( RQCORR*SGNDEF.GE.ZERO )
     $                       .AND.( LAMBDA + RQCORR.LE. RIGHT)
     $                       .AND.( LAMBDA + RQCORR.GE. LEFT)
     $                       ) THEN
                           USEDRQ = .TRUE.
<span class="comment">*</span><span class="comment">                          Store new midpoint of bisection interval in WORK
</span>                           IF(SGNDEF.EQ.ONE) THEN
<span class="comment">*</span><span class="comment">                             The current LAMBDA is on the left of the true
</span><span class="comment">*</span><span class="comment">                             eigenvalue
</span>                              LEFT = LAMBDA
<span class="comment">*</span><span class="comment">                             We prefer to assume that the error estimate
</span><span class="comment">*</span><span class="comment">                             is correct. We could make the interval not
</span><span class="comment">*</span><span class="comment">                             as a bracket but to be modified if the RQCORR
</span><span class="comment">*</span><span class="comment">                             chooses to. In this case, the RIGHT side should
</span><span class="comment">*</span><span class="comment">                             be modified as follows:
</span><span class="comment">*</span><span class="comment">                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
</span>                           ELSE
<span class="comment">*</span><span class="comment">                             The current LAMBDA is on the right of the true
</span><span class="comment">*</span><span class="comment">                             eigenvalue
</span>                              RIGHT = LAMBDA
<span class="comment">*</span><span class="comment">                             See comment about assuming the error estimate is
</span><span class="comment">*</span><span class="comment">                             correct above.
</span><span class="comment">*</span><span class="comment">                              LEFT = MIN(LEFT, LAMBDA + RQCORR)
</span>                           ENDIF
                           WORK( WINDEX ) =
     $                       HALF * (RIGHT + LEFT)
<span class="comment">*</span><span class="comment">                          Take RQCORR since it has the correct sign and
</span><span class="comment">*</span><span class="comment">                          improves the iterate reasonably
</span>                           LAMBDA = LAMBDA + RQCORR
<span class="comment">*</span><span class="comment">                          Update width of error interval
</span>                           WERR( WINDEX ) =
     $                             HALF * (RIGHT-LEFT)
                        ELSE
                           NEEDBS = .TRUE.
                        ENDIF
                        IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
<span class="comment">*</span><span class="comment">                             The eigenvalue is computed to bisection accuracy
</span><span class="comment">*</span><span class="comment">                             compute eigenvector and stop
</span>                           USEDBS = .TRUE.
                           GOTO 120
                        ELSEIF( ITER.LT.MAXITR ) THEN
                           GOTO 120
                        ELSEIF( ITER.EQ.MAXITR ) THEN
                           NEEDBS = .TRUE.
                           GOTO 120
                        ELSE
                           INFO = 5
                           RETURN
                        END IF
                     ELSE
                        STP2II = .FALSE.
        IF(USEDRQ .AND. USEDBS .AND.
     $                     BSTRES.LE.RESID) THEN
                           LAMBDA = BSTW
                           STP2II = .TRUE.
                        ENDIF
                        IF (STP2II) THEN
<span class="comment">*</span><span class="comment">                          improve error angle by second step
</span>                           CALL <a name="ZLAR1V.840"></a><a href="zlar1v.f.html#ZLAR1V.1">ZLAR1V</a>( IN, 1, IN, LAMBDA,
     $                          D( IBEGIN ), L( IBEGIN ),
     $                          WORK(INDLD+IBEGIN-1),
     $                          WORK(INDLLD+IBEGIN-1),
     $                          PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
     $                          .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
     $                          IWORK( IINDR+WINDEX ),
     $                          ISUPPZ( 2*WINDEX-1 ),
     $                          NRMINV, RESID, RQCORR, WORK( INDWRK ) )
                        ENDIF
                        WORK( WINDEX ) = LAMBDA
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Compute FP-vector support w.r.t. whole matrix
</span><span class="comment">*</span><span class="comment">
</span>                     ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
                     ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
                     ZFROM = ISUPPZ( 2*WINDEX-1 )
                     ZTO = ISUPPZ( 2*WINDEX )
                     ISUPMN = ISUPMN + OLDIEN
                     ISUPMX = ISUPMX + OLDIEN
<span class="comment">*</span><span class="comment">                    Ensure vector is ok if support in the RQI has changed
</span>