chegs2.f.html

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

         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.114"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CHEGS2.114"></a><a href="chegs2.f.html#CHEGS2.1">CHEGS2</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ITYPE.EQ.1 ) THEN
         IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute inv(U')*A*inv(U)
</span><span class="comment">*</span><span class="comment">
</span>            DO 10 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update the upper triangle of A(k:n,k:n)
</span><span class="comment">*</span><span class="comment">
</span>               AKK = A( K, K )
               BKK = B( K, K )
               AKK = AKK / BKK**2
               A( K, K ) = AKK
               IF( K.LT.N ) THEN
                  CALL CSSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
                  CT = -HALF*AKK
                  CALL <a name="CLACGV.134"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-K, A( K, K+1 ), LDA )
                  CALL <a name="CLACGV.135"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-K, B( K, K+1 ), LDB )
                  CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
     $                        LDA )
                  CALL CHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA,
     $                        B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
                  CALL CAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
     $                        LDA )
                  CALL <a name="CLACGV.142"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-K, B( K, K+1 ), LDB )
                  CALL CTRSV( UPLO, <span class="string">'Conjugate transpose'</span>, <span class="string">'Non-unit'</span>,
     $                        N-K, B( K+1, K+1 ), LDB, A( K, K+1 ),
     $                        LDA )
                  CALL <a name="CLACGV.146"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-K, A( K, K+1 ), LDA )
               END IF
   10       CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute inv(L)*A*inv(L')
</span><span class="comment">*</span><span class="comment">
</span>            DO 20 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update the lower triangle of A(k:n,k:n)
</span><span class="comment">*</span><span class="comment">
</span>               AKK = A( K, K )
               BKK = B( K, K )
               AKK = AKK / BKK**2
               A( K, K ) = AKK
               IF( K.LT.N ) THEN
                  CALL CSSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
                  CT = -HALF*AKK
                  CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
                  CALL CHER2( UPLO, N-K, -CONE, A( K+1, K ), 1,
     $                        B( K+1, K ), 1, A( K+1, K+1 ), LDA )
                  CALL CAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
                  CALL CTRSV( UPLO, <span class="string">'No transpose'</span>, <span class="string">'Non-unit'</span>, N-K,
     $                        B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
               END IF
   20       CONTINUE
         END IF
      ELSE
         IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute U*A*U'
</span><span class="comment">*</span><span class="comment">
</span>            DO 30 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update the upper triangle of A(1:k,1:k)
</span><span class="comment">*</span><span class="comment">
</span>               AKK = A( K, K )
               BKK = B( K, K )
               CALL CTRMV( UPLO, <span class="string">'No transpose'</span>, <span class="string">'Non-unit'</span>, K-1, B,
     $                     LDB, A( 1, K ), 1 )
               CT = HALF*AKK
               CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
               CALL CHER2( UPLO, K-1, CONE, A( 1, K ), 1, B( 1, K ), 1,
     $                     A, LDA )
               CALL CAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
               CALL CSSCAL( K-1, BKK, A( 1, K ), 1 )
               A( K, K ) = AKK*BKK**2
   30       CONTINUE
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute L'*A*L
</span><span class="comment">*</span><span class="comment">
</span>            DO 40 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update the lower triangle of A(1:k,1:k)
</span><span class="comment">*</span><span class="comment">
</span>               AKK = A( K, K )
               BKK = B( K, K )
               CALL <a name="CLACGV.204"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( K-1, A( K, 1 ), LDA )
               CALL CTRMV( UPLO, <span class="string">'Conjugate transpose'</span>, <span class="string">'Non-unit'</span>, K-1,
     $                     B, LDB, A( K, 1 ), LDA )
               CT = HALF*AKK
               CALL <a name="CLACGV.208"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( K-1, B( K, 1 ), LDB )
               CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
               CALL CHER2( UPLO, K-1, CONE, A( K, 1 ), LDA, B( K, 1 ),
     $                     LDB, A, LDA )
               CALL CAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
               CALL <a name="CLACGV.213"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( K-1, B( K, 1 ), LDB )
               CALL CSSCAL( K-1, BKK, A( K, 1 ), LDA )
               CALL <a name="CLACGV.215"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( K-1, A( K, 1 ), LDA )
               A( K, K ) = AKK*BKK**2
   40       CONTINUE
         END IF
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CHEGS2.222"></a><a href="chegs2.f.html#CHEGS2.1">CHEGS2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>