dlarre.f.html

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

</span><span class="comment">*</span><span class="comment">        Find local outer bounds GL,GU for the block
</span>         GL = D(IBEGIN)
         GU = D(IBEGIN)
         DO 15 I = IBEGIN , IEND
            GL = MIN( GERS( 2*I-1 ), GL )
            GU = MAX( GERS( 2*I ), GU )
 15      CONTINUE
         SPDIAM = GU - GL

         IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
<span class="comment">*</span><span class="comment">           Count the number of eigenvalues in the current block.
</span>            MB = 0
            DO 20 I = WBEGIN,MM
               IF( IBLOCK(I).EQ.JBLK ) THEN
                  MB = MB+1
               ELSE
                  GOTO 21
               ENDIF
 20         CONTINUE
 21         CONTINUE

            IF( MB.EQ.0) THEN
<span class="comment">*</span><span class="comment">              No eigenvalue in the current block lies in the desired range
</span><span class="comment">*</span><span class="comment">              E( IEND ) holds the shift for the initial RRR
</span>               E( IEND ) = ZERO
               IBEGIN = IEND + 1
               GO TO 170
            ELSE

<span class="comment">*</span><span class="comment">              Decide whether dqds or bisection is more efficient
</span>               USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
               WEND = WBEGIN + MB - 1
<span class="comment">*</span><span class="comment">              Calculate gaps for the current block
</span><span class="comment">*</span><span class="comment">              In later stages, when representations for individual
</span><span class="comment">*</span><span class="comment">              eigenvalues are different, we use SIGMA = E( IEND ).
</span>               SIGMA = ZERO
               DO 30 I = WBEGIN, WEND - 1
                  WGAP( I ) = MAX( ZERO,
     $                        W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
 30            CONTINUE
               WGAP( WEND ) = MAX( ZERO,
     $                     VU - SIGMA - (W( WEND )+WERR( WEND )))
<span class="comment">*</span><span class="comment">              Find local index of the first and last desired evalue.
</span>               INDL = INDEXW(WBEGIN)
               INDU = INDEXW( WEND )
            ENDIF
         ENDIF
         IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
<span class="comment">*</span><span class="comment">           Case of DQDS
</span><span class="comment">*</span><span class="comment">           Find approximations to the extremal eigenvalues of the block
</span>            CALL <a name="DLARRK.420"></a><a href="dlarrk.f.html#DLARRK.1">DLARRK</a>( IN, 1, GL, GU, D(IBEGIN),
     $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
            IF( IINFO.NE.0 ) THEN
               INFO = -1
               RETURN
            ENDIF
            ISLEFT = MAX(GL, TMP - TMP1
     $               - HNDRD * EPS* ABS(TMP - TMP1))

            CALL <a name="DLARRK.429"></a><a href="dlarrk.f.html#DLARRK.1">DLARRK</a>( IN, IN, GL, GU, D(IBEGIN),
     $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
            IF( IINFO.NE.0 ) THEN
               INFO = -1
               RETURN
            ENDIF
            ISRGHT = MIN(GU, TMP + TMP1
     $                 + HNDRD * EPS * ABS(TMP + TMP1))
<span class="comment">*</span><span class="comment">           Improve the estimate of the spectral diameter
</span>            SPDIAM = ISRGHT - ISLEFT
         ELSE
<span class="comment">*</span><span class="comment">           Case of bisection
</span><span class="comment">*</span><span class="comment">           Find approximations to the wanted extremal eigenvalues
</span>            ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
     $                  - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
            ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
     $                  + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
         ENDIF


<span class="comment">*</span><span class="comment">        Decide whether the base representation for the current block
</span><span class="comment">*</span><span class="comment">        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
</span><span class="comment">*</span><span class="comment">        should be on the left or the right end of the current block.
</span><span class="comment">*</span><span class="comment">        The strategy is to shift to the end which is &quot;more populated&quot;
</span><span class="comment">*</span><span class="comment">        Furthermore, decide whether to use DQDS for the computation of
</span><span class="comment">*</span><span class="comment">        the eigenvalue approximations at the end of <a name="DLARRE.454"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> or bisection.
</span><span class="comment">*</span><span class="comment">        dqds is chosen if all eigenvalues are desired or the number of
</span><span class="comment">*</span><span class="comment">        eigenvalues to be computed is large compared to the blocksize.
</span>         IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
<span class="comment">*</span><span class="comment">           If all the eigenvalues have to be computed, we use dqd
</span>            USEDQD = .TRUE.
<span class="comment">*</span><span class="comment">           INDL is the local index of the first eigenvalue to compute
</span>            INDL = 1
            INDU = IN
<span class="comment">*</span><span class="comment">           MB =  number of eigenvalues to compute
</span>            MB = IN
            WEND = WBEGIN + MB - 1
<span class="comment">*</span><span class="comment">           Define 1/4 and 3/4 points of the spectrum
</span>            S1 = ISLEFT + FOURTH * SPDIAM
            S2 = ISRGHT - FOURTH * SPDIAM
         ELSE
<span class="comment">*</span><span class="comment">           <a name="DLARRD.470"></a><a href="dlarrd.f.html#DLARRD.1">DLARRD</a> has computed IBLOCK and INDEXW for each eigenvalue
</span><span class="comment">*</span><span class="comment">           approximation.
</span><span class="comment">*</span><span class="comment">           choose sigma
</span>            IF( USEDQD ) THEN
               S1 = ISLEFT + FOURTH * SPDIAM
               S2 = ISRGHT - FOURTH * SPDIAM
            ELSE
               TMP = MIN(ISRGHT,VU) -  MAX(ISLEFT,VL)
               S1 =  MAX(ISLEFT,VL) + FOURTH * TMP
               S2 =  MIN(ISRGHT,VU) - FOURTH * TMP
            ENDIF
         ENDIF

<span class="comment">*</span><span class="comment">        Compute the negcount at the 1/4 and 3/4 points
</span>         IF(MB.GT.1) THEN
            CALL <a name="DLARRC.485"></a><a href="dlarrc.f.html#DLARRC.1">DLARRC</a>( <span class="string">'T'</span>, IN, S1, S2, D(IBEGIN),
     $                    E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
         ENDIF

         IF(MB.EQ.1) THEN
            SIGMA = GL
            SGNDEF = ONE
         ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
            IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
               SIGMA = MAX(ISLEFT,GL)
            ELSEIF( USEDQD ) THEN
<span class="comment">*</span><span class="comment">              use Gerschgorin bound as shift to get pos def matrix
</span><span class="comment">*</span><span class="comment">              for dqds
</span>               SIGMA = ISLEFT
            ELSE
<span class="comment">*</span><span class="comment">              use approximation of the first desired eigenvalue of the
</span><span class="comment">*</span><span class="comment">              block as shift
</span>               SIGMA = MAX(ISLEFT,VL)
            ENDIF
            SGNDEF = ONE
         ELSE
            IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
               SIGMA = MIN(ISRGHT,GU)
            ELSEIF( USEDQD ) THEN
<span class="comment">*</span><span class="comment">              use Gerschgorin bound as shift to get neg def matrix
</span><span class="comment">*</span><span class="comment">              for dqds
</span>               SIGMA = ISRGHT
            ELSE
<span class="comment">*</span><span class="comment">              use approximation of the first desired eigenvalue of the
</span><span class="comment">*</span><span class="comment">              block as shift
</span>               SIGMA = MIN(ISRGHT,VU)
            ENDIF
            SGNDEF = -ONE
         ENDIF


<span class="comment">*</span><span class="comment">        An initial SIGMA has been chosen that will be used for computing
</span><span class="comment">*</span><span class="comment">        T - SIGMA I = L D L^T
</span><span class="comment">*</span><span class="comment">        Define the increment TAU of the shift in case the initial shift
</span><span class="comment">*</span><span class="comment">        needs to be refined to obtain a factorization with not too much
</span><span class="comment">*</span><span class="comment">        element growth.
</span>         IF( USEDQD ) THEN
<span class="comment">*</span><span class="comment">           The initial SIGMA was to the outer end of the spectrum
</span><span class="comment">*</span><span class="comment">           the matrix is definite and we need not retreat.
</span>            TAU = SPDIAM*EPS*N + TWO*PIVMIN
         ELSE
            IF(MB.GT.1) THEN
               CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
               AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN))
               IF( SGNDEF.EQ.ONE ) THEN
                  TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
                  TAU = MAX(TAU,WERR(WBEGIN))
               ELSE
                  TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
                  TAU = MAX(TAU,WERR(WEND))
               ENDIF
            ELSE
               TAU = WERR(WBEGIN)
            ENDIF
         ENDIF
<span class="comment">*</span><span class="comment">
</span>         DO 80 IDUM = 1, MAXTRY
<span class="comment">*</span><span class="comment">           Compute L D L^T factorization of tridiagonal matrix T - sigma I.
</span><span class="comment">*</span><span class="comment">           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
</span><span class="comment">*</span><span class="comment">           pivots in WORK(2*IN+1:3*IN)
</span>            DPIVOT = D( IBEGIN ) - SIGMA
            WORK( 1 ) = DPIVOT
            DMAX = ABS( WORK(1) )
            J = IBEGIN
            DO 70 I = 1, IN - 1
               WORK( 2*IN+I ) = ONE / WORK( I )
               TMP = E( J )*WORK( 2*IN+I )
               WORK( IN+I ) = TMP
               DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
               WORK( I+1 ) = DPIVOT
               DMAX = MAX( DMAX, ABS(DPIVOT) )
               J = J + 1
 70         CONTINUE
<span class="comment">*</span><span class="comment">           check for element growth