zstemr.f.html

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

</span>         IINFO = -1
      ENDIF
<span class="comment">*</span><span class="comment">     Set the splitting criterion
</span>      IF (IINFO.EQ.0) THEN
         THRESH = EPS
      ELSE
         THRESH = -EPS
<span class="comment">*</span><span class="comment">        relative accuracy is desired but T does not guarantee it
</span>         TRYRAC = .FALSE.
      ENDIF
<span class="comment">*</span><span class="comment">
</span>      IF( TRYRAC ) THEN
<span class="comment">*</span><span class="comment">        Copy original diagonal, needed to guarantee relative accuracy
</span>         CALL DCOPY(N,D,1,WORK(INDD),1)
      ENDIF
<span class="comment">*</span><span class="comment">     Store the squares of the offdiagonal values of T
</span>      DO 5 J = 1, N-1
         WORK( INDE2+J-1 ) = E(J)**2
 5    CONTINUE

<span class="comment">*</span><span class="comment">     Set the tolerance parameters for bisection
</span>      IF( .NOT.WANTZ ) THEN
<span class="comment">*</span><span class="comment">        <a name="DLARRE.523"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> computes the eigenvalues to full precision.
</span>         RTOL1 = FOUR * EPS
         RTOL2 = FOUR * EPS
      ELSE
<span class="comment">*</span><span class="comment">        <a name="DLARRE.527"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> computes the eigenvalues to less than full precision.
</span><span class="comment">*</span><span class="comment">        <a name="ZLARRV.528"></a><a href="zlarrv.f.html#ZLARRV.1">ZLARRV</a> will refine the eigenvalue approximations, and we only
</span><span class="comment">*</span><span class="comment">        need less accurate initial bisection in <a name="DLARRE.529"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>.
</span><span class="comment">*</span><span class="comment">        Note: these settings do only affect the subset case and <a name="DLARRE.530"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>
</span>         RTOL1 = SQRT(EPS)
         RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
      ENDIF
      CALL <a name="DLARRE.534"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a>( RANGE, N, WL, WU, IIL, IIU, D, E,
     $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
     $             IWORK( IINSPL ), M, W, WORK( INDERR ),
     $             WORK( INDGP ), IWORK( IINDBL ),
     $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
     $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )
      IF( IINFO.NE.0 ) THEN
         INFO = 10 + ABS( IINFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">     Note that if RANGE .NE. 'V', <a name="DLARRE.544"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> computes bounds on the desired
</span><span class="comment">*</span><span class="comment">     part of the spectrum. All desired eigenvalues are contained in
</span><span class="comment">*</span><span class="comment">     (WL,WU]
</span>

      IF( WANTZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the desired eigenvectors corresponding to the computed
</span><span class="comment">*</span><span class="comment">        eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLARRV.554"></a><a href="zlarrv.f.html#ZLARRV.1">ZLARRV</a>( N, WL, WU, D, E,
     $                PIVMIN, IWORK( IINSPL ), M,
     $                1, M, MINRGP, RTOL1, RTOL2,
     $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
     $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
     $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
         IF( IINFO.NE.0 ) THEN
            INFO = 20 + ABS( IINFO )
            RETURN
         END IF
      ELSE
<span class="comment">*</span><span class="comment">        <a name="DLARRE.565"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> computes eigenvalues of the (shifted) root representation
</span><span class="comment">*</span><span class="comment">        <a name="ZLARRV.566"></a><a href="zlarrv.f.html#ZLARRV.1">ZLARRV</a> returns the eigenvalues of the unshifted matrix.
</span><span class="comment">*</span><span class="comment">        However, if the eigenvectors are not desired by the user, we need
</span><span class="comment">*</span><span class="comment">        to apply the corresponding shifts from <a name="DLARRE.568"></a><a href="dlarre.f.html#DLARRE.1">DLARRE</a> to obtain the
</span><span class="comment">*</span><span class="comment">        eigenvalues of the original matrix.
</span>         DO 20 J = 1, M
            ITMP = IWORK( IINDBL+J-1 )
            W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
 20      CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>
      IF ( TRYRAC ) THEN
<span class="comment">*</span><span class="comment">        Refine computed eigenvalues so that they are relatively accurate
</span><span class="comment">*</span><span class="comment">        with respect to the original matrix T.
</span>         IBEGIN = 1
         WBEGIN = 1
         DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
            IEND = IWORK( IINSPL+JBLK-1 )
            IN = IEND - IBEGIN + 1
            WEND = WBEGIN - 1
<span class="comment">*</span><span class="comment">           check if any eigenvalues have to be refined in this block
</span> 36         CONTINUE
            IF( WEND.LT.M ) THEN
               IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
                  WEND = WEND + 1
                  GO TO 36
               END IF
            END IF
            IF( WEND.LT.WBEGIN ) THEN
               IBEGIN = IEND + 1
               GO TO 39
            END IF

            OFFSET = IWORK(IINDW+WBEGIN-1)-1
            IFIRST = IWORK(IINDW+WBEGIN-1)
            ILAST = IWORK(IINDW+WEND-1)
            RTOL2 = FOUR * EPS
            CALL <a name="DLARRJ.603"></a><a href="dlarrj.f.html#DLARRJ.1">DLARRJ</a>( IN,
     $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
     $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
     $                   WORK( INDERR+WBEGIN-1 ),
     $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
     $                   TNRM, IINFO )
            IBEGIN = IEND + 1
            WBEGIN = WEND + 1
 39      CONTINUE
      ENDIF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If matrix was scaled, then rescale eigenvalues appropriately.
</span><span class="comment">*</span><span class="comment">
</span>      IF( SCALE.NE.ONE ) THEN
         CALL DSCAL( M, ONE / SCALE, W, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If eigenvalues are not in increasing order, then sort them,
</span><span class="comment">*</span><span class="comment">     possibly along with eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>      IF( NSPLIT.GT.1 ) THEN
         IF( .NOT. WANTZ ) THEN
            CALL <a name="DLASRT.625"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>( <span class="string">'I'</span>, M, W, IINFO )
            IF( IINFO.NE.0 ) THEN
               INFO = 3
               RETURN
            END IF
         ELSE
            DO 60 J = 1, M - 1
               I = 0
               TMP = W( J )
               DO 50 JJ = J + 1, M
                  IF( W( JJ ).LT.TMP ) THEN
                     I = JJ
                     TMP = W( JJ )
                  END IF
 50            CONTINUE
               IF( I.NE.0 ) THEN
                  W( I ) = W( J )
                  W( J ) = TMP
                  IF( WANTZ ) THEN
                     CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
                     ITMP = ISUPPZ( 2*I-1 )
                     ISUPPZ( 2*I-1 ) = ISUPPZ( 2*J-1 )
                     ISUPPZ( 2*J-1 ) = ITMP
                     ITMP = ISUPPZ( 2*I )
                     ISUPPZ( 2*I ) = ISUPPZ( 2*J )
                     ISUPPZ( 2*J ) = ITMP
                  END IF
               END IF
 60         CONTINUE
         END IF
      ENDIF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWMIN
      IWORK( 1 ) = LIWMIN
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZSTEMR.662"></a><a href="zstemr.f.html#ZSTEMR.1">ZSTEMR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>