chegvd.f.html

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

</span><span class="comment">*</span><span class="comment">                    minor of order i of B is not positive definite.
</span><span class="comment">*</span><span class="comment">                    The factorization of B could not be completed and
</span><span class="comment">*</span><span class="comment">                    no eigenvalues or eigenvectors were computed.
</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">     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Modified so that no backsubstitution is performed if <a name="CHEEVD.158"></a><a href="cheevd.f.html#CHEEVD.1">CHEEVD</a> fails to
</span><span class="comment">*</span><span class="comment">  converge (NEIG in old code could be greater than N causing out of
</span><span class="comment">*</span><span class="comment">  bounds reference to A - reported by Ralf Meyer).  Also corrected the
</span><span class="comment">*</span><span class="comment">  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
</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>      COMPLEX            CONE
      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, UPPER, WANTZ
      CHARACTER          TRANS
      INTEGER            LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.174"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      EXTERNAL           <a name="LSAME.175"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="CHEEVD.178"></a><a href="cheevd.f.html#CHEEVD.1">CHEEVD</a>, <a name="CHEGST.178"></a><a href="chegst.f.html#CHEGST.1">CHEGST</a>, <a name="CPOTRF.178"></a><a href="cpotrf.f.html#CPOTRF.1">CPOTRF</a>, CTRMM, CTRSM, <a name="XERBLA.178"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX, REAL
<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">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      WANTZ = <a name="LSAME.187"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBZ, <span class="string">'V'</span> )
      UPPER = <a name="LSAME.188"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( N.LE.1 ) THEN
         LWMIN = 1
         LRWMIN = 1
         LIWMIN = 1
      ELSE IF( WANTZ ) THEN
         LWMIN = 2*N + N*N
         LRWMIN = 1 + 5*N + 2*N*N
         LIWMIN = 3 + 5*N
      ELSE
         LWMIN = N + 1
         LRWMIN = N
         LIWMIN = 1
      END IF
      LOPT = LWMIN
      LROPT = LRWMIN
      LIOPT = LIWMIN
      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
         INFO = -1
      ELSE IF( .NOT.( WANTZ .OR. <a name="LSAME.210"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBZ, <span class="string">'N'</span> ) ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( UPPER .OR. <a name="LSAME.212"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) ) ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -8
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = LOPT
         RWORK( 1 ) = LROPT
         IWORK( 1 ) = LIOPT
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -11
         ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -13
         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
            INFO = -15
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.237"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CHEGVD.237"></a><a href="chegvd.f.html#CHEGVD.1">CHEGVD</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Form a Cholesky factorization of B.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CPOTRF.250"></a><a href="cpotrf.f.html#CPOTRF.1">CPOTRF</a>( UPLO, N, B, LDB, INFO )
      IF( INFO.NE.0 ) THEN
         INFO = N + INFO
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Transform problem to standard eigenvalue problem and solve.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CHEGST.258"></a><a href="chegst.f.html#CHEGST.1">CHEGST</a>( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
      CALL <a name="CHEEVD.259"></a><a href="cheevd.f.html#CHEEVD.1">CHEEVD</a>( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK,
     $             IWORK, LIWORK, INFO )
      LOPT = MAX( REAL( LOPT ), REAL( WORK( 1 ) ) )
      LROPT = MAX( REAL( LROPT ), REAL( RWORK( 1 ) ) )
      LIOPT = MAX( REAL( LIOPT ), REAL( IWORK( 1 ) ) )
<span class="comment">*</span><span class="comment">
</span>      IF( WANTZ .AND. INFO.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Backtransform eigenvectors to the original problem.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
</span><span class="comment">*</span><span class="comment">           backtransform eigenvectors: x = inv(L)'*y or inv(U)*y
</span><span class="comment">*</span><span class="comment">
</span>            IF( UPPER ) THEN
               TRANS = <span class="string">'N'</span>
            ELSE
               TRANS = <span class="string">'C'</span>
            END IF
<span class="comment">*</span><span class="comment">
</span>            CALL CTRSM( <span class="string">'Left'</span>, UPLO, TRANS, <span class="string">'Non-unit'</span>, N, N, CONE,
     $                  B, LDB, A, LDA )
<span class="comment">*</span><span class="comment">
</span>         ELSE IF( ITYPE.EQ.3 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           For B*A*x=(lambda)*x;
</span><span class="comment">*</span><span class="comment">           backtransform eigenvectors: x = L*y or U'*y
</span><span class="comment">*</span><span class="comment">
</span>            IF( UPPER ) THEN
               TRANS = <span class="string">'C'</span>
            ELSE
               TRANS = <span class="string">'N'</span>
            END IF
<span class="comment">*</span><span class="comment">
</span>            CALL CTRMM( <span class="string">'Left'</span>, UPLO, TRANS, <span class="string">'Non-unit'</span>, N, N, CONE,
     $                  B, LDB, A, LDA )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LOPT
      RWORK( 1 ) = LROPT
      IWORK( 1 ) = LIOPT
<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="CHEGVD.305"></a><a href="chegvd.f.html#CHEGVD.1">CHEGVD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>