dsyevr.f.html

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

</span>      INDISP = INDIBL + N
<span class="comment">*</span><span class="comment">     IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
</span><span class="comment">*</span><span class="comment">     that corresponding to eigenvectors that fail to converge in
</span><span class="comment">*</span><span class="comment">     <a name="DSTEIN.417"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.  This information is discarded; if any fail, the driver
</span><span class="comment">*</span><span class="comment">     returns INFO &gt; 0.
</span>      INDIFL = INDISP + N
<span class="comment">*</span><span class="comment">     INDIWO is the offset of the remaining integer workspace.
</span>      INDIWO = INDISP + N

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Call <a name="DSYTRD.424"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a> to reduce symmetric matrix to tridiagonal form.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DSYTRD.426"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>( UPLO, N, A, LDA, WORK( INDD ), WORK( INDE ),
     $             WORK( INDTAU ), WORK( INDWK ), LLWORK, IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If all eigenvalues are desired
</span><span class="comment">*</span><span class="comment">     then call <a name="DSTERF.430"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> or <a name="DSTEMR.430"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> and <a name="DORMTR.430"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( ALLEIG .OR. ( INDEIG .AND. IL.EQ.1 .AND. IU.EQ.N ) ) .AND.
     $    IEEEOK.EQ.1 ) THEN
         IF( .NOT.WANTZ ) THEN
            CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
            CALL <a name="DSTERF.437"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>( N, W, WORK( INDEE ), INFO )
         ELSE
            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
            CALL DCOPY( N, WORK( INDD ), 1, WORK( INDDD ), 1 )
<span class="comment">*</span><span class="comment">
</span>            IF (ABSTOL .LE. TWO*N*EPS) THEN
               TRYRAC = .TRUE.
            ELSE
               TRYRAC = .FALSE.
            END IF
            CALL <a name="DSTEMR.447"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>( JOBZ, <span class="string">'A'</span>, N, WORK( INDDD ), WORK( INDEE ),
     $                   VL, VU, IL, IU, M, W, Z, LDZ, N, ISUPPZ,
     $                   TRYRAC, WORK( INDWK ), LWORK, IWORK, LIWORK,
     $                   INFO )
<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">        Apply orthogonal matrix used in reduction to tridiagonal
</span><span class="comment">*</span><span class="comment">        form to eigenvectors returned by <a name="DSTEIN.455"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.
</span><span class="comment">*</span><span class="comment">
</span>            IF( WANTZ .AND. INFO.EQ.0 ) THEN
               INDWKN = INDE
               LLWRKN = LWORK - INDWKN + 1
               CALL <a name="DORMTR.460"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>( <span class="string">'L'</span>, UPLO, <span class="string">'N'</span>, N, M, A, LDA,
     $                      WORK( INDTAU ), Z, LDZ, WORK( INDWKN ),
     $                      LLWRKN, IINFO )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>         IF( INFO.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">           Everything worked.  Skip <a name="DSTEBZ.468"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>/<a name="DSTEIN.468"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.  IWORK(:) are
</span><span class="comment">*</span><span class="comment">           undefined.
</span>            M = N
            GO TO 30
         END IF
         INFO = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Otherwise, call <a name="DSTEBZ.476"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and, if eigenvectors are desired, <a name="DSTEIN.476"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.
</span><span class="comment">*</span><span class="comment">     Also call <a name="DSTEBZ.477"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and <a name="DSTEIN.477"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a> if <a name="DSTEMR.477"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a> fails.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         ORDER = <span class="string">'B'</span>
      ELSE
         ORDER = <span class="string">'E'</span>
      END IF

      CALL <a name="DSTEBZ.485"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
     $             WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
     $             IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWK ),
     $             IWORK( INDIWO ), INFO )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         CALL <a name="DSTEIN.491"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>( N, WORK( INDD ), WORK( INDE ), M, W,
     $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
     $                WORK( INDWK ), IWORK( INDIWO ), IWORK( INDIFL ),
     $                INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply orthogonal matrix used in reduction to tridiagonal
</span><span class="comment">*</span><span class="comment">        form to eigenvectors returned by <a name="DSTEIN.497"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>.
</span><span class="comment">*</span><span class="comment">
</span>         INDWKN = INDE
         LLWRKN = LWORK - INDWKN + 1
         CALL <a name="DORMTR.501"></a><a href="dormtr.f.html#DORMTR.1">DORMTR</a>( <span class="string">'L'</span>, UPLO, <span class="string">'N'</span>, N, M, A, LDA, WORK( INDTAU ), Z,
     $                LDZ, WORK( INDWKN ), LLWRKN, IINFO )
      END IF
<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><span class="comment">*</span><span class="comment">  Jump here if <a name="DSTEMR.507"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>/<a name="DSTEIN.507"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a> succeeded.
</span>   30 CONTINUE
      IF( ISCALE.EQ.1 ) THEN
         IF( INFO.EQ.0 ) THEN
            IMAX = M
         ELSE
            IMAX = INFO - 1
         END IF
         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If eigenvalues are not in order, then sort them, along with
</span><span class="comment">*</span><span class="comment">     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK.
</span><span class="comment">*</span><span class="comment">     It may not be initialized (if <a name="DSTEMR.520"></a><a href="dstemr.f.html#DSTEMR.1">DSTEMR</a>/<a name="DSTEIN.520"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a> succeeded), and we do
</span><span class="comment">*</span><span class="comment">     not return this detailed information to the user.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         DO 50 J = 1, M - 1
            I = 0
            TMP1 = W( J )
            DO 40 JJ = J + 1, M
               IF( W( JJ ).LT.TMP1 ) THEN
                  I = JJ
                  TMP1 = W( JJ )
               END IF
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( I.NE.0 ) THEN
               W( I ) = W( J )
               W( J ) = TMP1
               CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
            END IF
   50    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set WORK(1) to optimal workspace size.
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWKOPT
      IWORK( 1 ) = LIWMIN
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DSYEVR.549"></a><a href="dsyevr.f.html#DSYEVR.1">DSYEVR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>