dsyevx.f.html

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

         SIGMA = RMAX / ANRM
      END IF
      IF( ISCALE.EQ.1 ) THEN
         IF( LOWER ) THEN
            DO 10 J = 1, N
               CALL DSCAL( N-J+1, SIGMA, A( J, J ), 1 )
   10       CONTINUE
         ELSE
            DO 20 J = 1, N
               CALL DSCAL( J, SIGMA, A( 1, J ), 1 )
   20       CONTINUE
         END IF
         IF( ABSTOL.GT.0 )
     $      ABSTLL = ABSTOL*SIGMA
         IF( VALEIG ) THEN
            VLL = VL*SIGMA
            VUU = VU*SIGMA
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Call <a name="DSYTRD.308"></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>      INDTAU = 1
      INDE = INDTAU + N
      INDD = INDE + N
      INDWRK = INDD + N
      LLWORK = LWORK - INDWRK + 1
      CALL <a name="DSYTRD.315"></a><a href="dsytrd.f.html#DSYTRD.1">DSYTRD</a>( UPLO, N, A, LDA, WORK( INDD ), WORK( INDE ),
     $             WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If all eigenvalues are desired and ABSTOL is less than or equal to
</span><span class="comment">*</span><span class="comment">     zero, then call <a name="DSTERF.319"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> or <a name="DORGTR.319"></a><a href="dorgtr.f.html#DORGTR.1">DORGTR</a> and <a name="SSTEQR.319"></a><a href="ssteqr.f.html#SSTEQR.1">SSTEQR</a>.  If this fails for
</span><span class="comment">*</span><span class="comment">     some eigenvalue, then try <a name="DSTEBZ.320"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a>.
</span><span class="comment">*</span><span class="comment">
</span>      TEST = .FALSE.
      IF( INDEIG ) THEN
         IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
            TEST = .TRUE.
         END IF
      END IF
      IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
         CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
         INDEE = INDWRK + 2*N
         IF( .NOT.WANTZ ) THEN
            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
            CALL <a name="DSTERF.333"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>( N, W, WORK( INDEE ), INFO )
         ELSE
            CALL <a name="DLACPY.335"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'A'</span>, N, N, A, LDA, Z, LDZ )
            CALL <a name="DORGTR.336"></a><a href="dorgtr.f.html#DORGTR.1">DORGTR</a>( UPLO, N, Z, LDZ, WORK( INDTAU ),
     $                   WORK( INDWRK ), LLWORK, IINFO )
            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
            CALL <a name="DSTEQR.339"></a><a href="dsteqr.f.html#DSTEQR.1">DSTEQR</a>( JOBZ, N, W, WORK( INDEE ), Z, LDZ,
     $                   WORK( INDWRK ), INFO )
            IF( INFO.EQ.0 ) THEN
               DO 30 I = 1, N
                  IFAIL( I ) = 0
   30          CONTINUE
            END IF
         END IF
         IF( INFO.EQ.0 ) THEN
            M = N
            GO TO 40
         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.354"></a><a href="dstebz.f.html#DSTEBZ.1">DSTEBZ</a> and, if eigenvectors are desired, <a name="SSTEIN.354"></a><a href="sstein.f.html#SSTEIN.1">SSTEIN</a>.
</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
      INDIBL = 1
      INDISP = INDIBL + N
      INDIWO = INDISP + N
      CALL <a name="DSTEBZ.364"></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( INDWRK ),
     $             IWORK( INDIWO ), INFO )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         CALL <a name="DSTEIN.370"></a><a href="dstein.f.html#DSTEIN.1">DSTEIN</a>( N, WORK( INDD ), WORK( INDE ), M, W,
     $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
     $                WORK( INDWRK ), IWORK( INDIWO ), IFAIL, 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.375"></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.379"></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>   40 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.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTZ ) THEN
         DO 60 J = 1, M - 1
            I = 0
            TMP1 = W( J )
            DO 50 JJ = J + 1, M
               IF( W( JJ ).LT.TMP1 ) THEN
                  I = JJ
                  TMP1 = W( JJ )
               END IF
   50       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( I.NE.0 ) THEN
               ITMP1 = IWORK( INDIBL+I-1 )
               W( I ) = W( J )
               IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
               W( J ) = TMP1
               IWORK( INDIBL+J-1 ) = ITMP1
               CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
               IF( INFO.NE.0 ) THEN
                  ITMP1 = IFAIL( I )
                  IFAIL( I ) = IFAIL( J )
                  IFAIL( J ) = ITMP1
               END IF
            END IF
   60    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
<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="DSYEVX.431"></a><a href="dsyevx.f.html#DSYEVX.1">DSYEVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>