zlarrv.f.html

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

</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Beresford Parlett, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Jim Demmel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Inderjit Dhillon, University of Texas, Austin, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">     Christof Voemel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      INTEGER            MAXITR
      PARAMETER          ( MAXITR = 10 )
      COMPLEX*16         CZERO
      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
      DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, HALF
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
     $                     TWO = 2.0D0, THREE = 3.0D0,
     $                     FOUR = 4.0D0, HALF = 0.5D0)
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
      INTEGER            DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
     $                   IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
     $                   INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
     $                   ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
     $                   NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
     $                   NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
     $                   OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
     $                   WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
     $                   ZUSEDW
      INTEGER            INDIN1, INDIN2
      DOUBLE PRECISION   BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
     $                   LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
     $                   RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
     $                   SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.200"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.201"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DCOPY, <a name="DLARRB.204"></a><a href="dlarrb.f.html#DLARRB.1">DLARRB</a>, <a name="DLARRF.204"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>, ZDSCAL, <a name="ZLAR1V.204"></a><a href="zlar1v.f.html#ZLAR1V.1">ZLAR1V</a>,
     $                   <a name="ZLASET.205"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC ABS, DBLE, MAX, MIN
      INTRINSIC DCMPLX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">     ..
</span>
<span class="comment">*</span><span class="comment">     The first N entries of WORK are reserved for the eigenvalues
</span>      INDLD = N+1
      INDLLD= 2*N+1
      INDIN1 = 3*N + 1
      INDIN2 = 4*N + 1
      INDWRK = 5*N + 1
      MINWSIZE = 12 * N

      DO 5 I= 1,MINWSIZE
         WORK( I ) = ZERO
 5    CONTINUE

<span class="comment">*</span><span class="comment">     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
</span><span class="comment">*</span><span class="comment">     factorization used to compute the FP vector
</span>      IINDR = 0
<span class="comment">*</span><span class="comment">     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
</span><span class="comment">*</span><span class="comment">     layer and the one above.
</span>      IINDC1 = N
      IINDC2 = 2*N
      IINDWK = 3*N + 1

      MINIWSIZE = 7 * N
      DO 10 I= 1,MINIWSIZE
         IWORK( I ) = 0
 10   CONTINUE

      ZUSEDL = 1
      IF(DOL.GT.1) THEN
<span class="comment">*</span><span class="comment">        Set lower bound for use of Z
</span>         ZUSEDL = DOL-1
      ENDIF
      ZUSEDU = M
      IF(DOU.LT.M) THEN
<span class="comment">*</span><span class="comment">        Set lower bound for use of Z
</span>         ZUSEDU = DOU+1
      ENDIF
<span class="comment">*</span><span class="comment">     The width of the part of Z that is used
</span>      ZUSEDW = ZUSEDU - ZUSEDL + 1


      CALL <a name="ZLASET.254"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, N, ZUSEDW, CZERO, CZERO,
     $                    Z(1,ZUSEDL), LDZ )

      EPS = <a name="DLAMCH.257"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      RQTOL = TWO * EPS
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set expert flags for standard code.
</span>      TRYRQC = .TRUE.

      IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
      ELSE
<span class="comment">*</span><span class="comment">        Only selected eigenpairs are computed. Since the other evalues
</span><span class="comment">*</span><span class="comment">        are not refined by RQ iteration, bisection has to compute to full
</span><span class="comment">*</span><span class="comment">        accuracy.
</span>         RTOL1 = FOUR * EPS
         RTOL2 = FOUR * EPS
      ENDIF

<span class="comment">*</span><span class="comment">     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
</span><span class="comment">*</span><span class="comment">     desired eigenvalues. The support of the nonzero eigenvector
</span><span class="comment">*</span><span class="comment">     entries is contained in the interval IBEGIN:IEND.
</span><span class="comment">*</span><span class="comment">     Remark that if k eigenpairs are desired, then the eigenvectors
</span><span class="comment">*</span><span class="comment">     are stored in k contiguous columns of Z.
</span>
<span class="comment">*</span><span class="comment">     DONE is the number of eigenvectors already computed
</span>      DONE = 0
      IBEGIN = 1
      WBEGIN = 1
      DO 170 JBLK = 1, IBLOCK( M )
         IEND = ISPLIT( JBLK )
         SIGMA = L( IEND )
<span class="comment">*</span><span class="comment">        Find the eigenvectors of the submatrix indexed IBEGIN
</span><span class="comment">*</span><span class="comment">        through IEND.
</span>         WEND = WBEGIN - 1
 15      CONTINUE
         IF( WEND.LT.M ) THEN
            IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
               WEND = WEND + 1
               GO TO 15
            END IF
         END IF
         IF( WEND.LT.WBEGIN ) THEN
            IBEGIN = IEND + 1
            GO TO 170
         ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
            IBEGIN = IEND + 1
            WBEGIN = WEND + 1
            GO TO 170
         END IF

<span class="comment">*</span><span class="comment">        Find local spectral diameter of the block
</span>         GL = GERS( 2*IBEGIN-1 )
         GU = GERS( 2*IBEGIN )
         DO 20 I = IBEGIN+1 , IEND
            GL = MIN( GERS( 2*I-1 ), GL )
            GU = MAX( GERS( 2*I ), GU )
 20      CONTINUE
         SPDIAM = GU - GL

<span class="comment">*</span><span class="comment">        OLDIEN is the last index of the previous block
</span>         OLDIEN = IBEGIN - 1
<span class="comment">*</span><span class="comment">        Calculate the size of the current block
</span>         IN = IEND - IBEGIN + 1
<span class="comment">*</span><span class="comment">        The number of eigenvalues in the current block
</span>         IM = WEND - WBEGIN + 1

<span class="comment">*</span><span class="comment">        This is for a 1x1 block
</span>         IF( IBEGIN.EQ.IEND ) THEN
            DONE = DONE+1
            Z( IBEGIN, WBEGIN ) = DCMPLX( ONE, ZERO )
            ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
            ISUPPZ( 2*WBEGIN ) = IBEGIN
            W( WBEGIN ) = W( WBEGIN ) + SIGMA
            WORK( WBEGIN ) = W( WBEGIN )
            IBEGIN = IEND + 1
            WBEGIN = WBEGIN + 1
            GO TO 170
         END IF

<span class="comment">*</span><span class="comment">        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
</span><span class="comment">*</span><span class="comment">        Note that these can be approximations, in this case, the corresp.
</span><span class="comment">*</span><span class="comment">        entries of WERR give the size of the uncertainty interval.