chetrs.f.html

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

HTML
418 行
字号
12
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.GT.N )
     $      GO TO 50
<span class="comment">*</span><span class="comment">
</span>         IF( IPIV( K ).GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           1 x 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by inv(U'(K)), where U(K) is the transformation
</span><span class="comment">*</span><span class="comment">           stored in column K of A.
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.GT.1 ) THEN
               CALL <a name="CLACGV.203"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL CGEMV( <span class="string">'Conjugate transpose'</span>, K-1, NRHS, -ONE, B,
     $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
               CALL <a name="CLACGV.206"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Interchange rows K and IPIV(K).
</span><span class="comment">*</span><span class="comment">
</span>            KP = IPIV( K )
            IF( KP.NE.K )
     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
            K = K + 1
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           2 x 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
</span><span class="comment">*</span><span class="comment">           stored in columns K and K+1 of A.
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.GT.1 ) THEN
               CALL <a name="CLACGV.223"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL CGEMV( <span class="string">'Conjugate transpose'</span>, K-1, NRHS, -ONE, B,
     $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
               CALL <a name="CLACGV.226"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLACGV.228"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K+1, 1 ), LDB )
               CALL CGEMV( <span class="string">'Conjugate transpose'</span>, K-1, NRHS, -ONE, B,
     $                     LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
               CALL <a name="CLACGV.231"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K+1, 1 ), LDB )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Interchange rows K and -IPIV(K).
</span><span class="comment">*</span><span class="comment">
</span>            KP = -IPIV( K )
            IF( KP.NE.K )
     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
            K = K + 2
         END IF
<span class="comment">*</span><span class="comment">
</span>         GO TO 40
   50    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve A*X = B, where A = L*D*L'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        First solve L*D*X = B, overwriting B with X.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        K is the main loop index, increasing from 1 to N in steps of
</span><span class="comment">*</span><span class="comment">        1 or 2, depending on the size of the diagonal blocks.
</span><span class="comment">*</span><span class="comment">
</span>         K = 1
   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If K &gt; N, exit from loop.
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.GT.N )
     $      GO TO 80
<span class="comment">*</span><span class="comment">
</span>         IF( IPIV( K ).GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           1 x 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Interchange rows K and IPIV(K).
</span><span class="comment">*</span><span class="comment">
</span>            KP = IPIV( K )
            IF( KP.NE.K )
     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by inv(L(K)), where L(K) is the transformation
</span><span class="comment">*</span><span class="comment">           stored in column K of A.
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.LT.N )
     $         CALL CGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
     $                     LDB, B( K+1, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by the inverse of the diagonal block.
</span><span class="comment">*</span><span class="comment">
</span>            S = REAL( ONE ) / REAL( A( K, K ) )
            CALL CSSCAL( NRHS, S, B( K, 1 ), LDB )
            K = K + 1
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           2 x 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Interchange rows K+1 and -IPIV(K).
</span><span class="comment">*</span><span class="comment">
</span>            KP = -IPIV( K )
            IF( KP.NE.K+1 )
     $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by inv(L(K)), where L(K) is the transformation
</span><span class="comment">*</span><span class="comment">           stored in columns K and K+1 of A.
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.LT.N-1 ) THEN
               CALL CGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
     $                     LDB, B( K+2, 1 ), LDB )
               CALL CGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
     $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by the inverse of the diagonal block.
</span><span class="comment">*</span><span class="comment">
</span>            AKM1K = A( K+1, K )
            AKM1 = A( K, K ) / CONJG( AKM1K )
            AK = A( K+1, K+1 ) / AKM1K
            DENOM = AKM1*AK - ONE
            DO 70 J = 1, NRHS
               BKM1 = B( K, J ) / CONJG( AKM1K )
               BK = B( K+1, J ) / AKM1K
               B( K, J ) = ( AK*BKM1-BK ) / DENOM
               B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
   70       CONTINUE
            K = K + 2
         END IF
<span class="comment">*</span><span class="comment">
</span>         GO TO 60
   80    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Next solve L'*X = B, overwriting B with X.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        K is the main loop index, decreasing from N to 1 in steps of
</span><span class="comment">*</span><span class="comment">        1 or 2, depending on the size of the diagonal blocks.
</span><span class="comment">*</span><span class="comment">
</span>         K = N
   90    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If K &lt; 1, exit from loop.
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.LT.1 )
     $      GO TO 100
<span class="comment">*</span><span class="comment">
</span>         IF( IPIV( K ).GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           1 x 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by inv(L'(K)), where L(K) is the transformation
</span><span class="comment">*</span><span class="comment">           stored in column K of A.
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.LT.N ) THEN
               CALL <a name="CLACGV.343"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL CGEMV( <span class="string">'Conjugate transpose'</span>, N-K, NRHS, -ONE,
     $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
     $                     B( K, 1 ), LDB )
               CALL <a name="CLACGV.347"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Interchange rows K and IPIV(K).
</span><span class="comment">*</span><span class="comment">
</span>            KP = IPIV( K )
            IF( KP.NE.K )
     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
            K = K - 1
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           2 x 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
</span><span class="comment">*</span><span class="comment">           stored in columns K-1 and K of A.
</span><span class="comment">*</span><span class="comment">
</span>            IF( K.LT.N ) THEN
               CALL <a name="CLACGV.364"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL CGEMV( <span class="string">'Conjugate transpose'</span>, N-K, NRHS, -ONE,
     $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
     $                     B( K, 1 ), LDB )
               CALL <a name="CLACGV.368"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLACGV.370"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K-1, 1 ), LDB )
               CALL CGEMV( <span class="string">'Conjugate transpose'</span>, N-K, NRHS, -ONE,
     $                     B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
     $                     B( K-1, 1 ), LDB )
               CALL <a name="CLACGV.374"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NRHS, B( K-1, 1 ), LDB )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Interchange rows K and -IPIV(K).
</span><span class="comment">*</span><span class="comment">
</span>            KP = -IPIV( K )
            IF( KP.NE.K )
     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
            K = K - 2
         END IF
<span class="comment">*</span><span class="comment">
</span>         GO TO 90
  100    CONTINUE
      END IF
<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="CHETRS.391"></a><a href="chetrs.f.html#CHETRS.1">CHETRS</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
12

chetrs.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 chetrs.f.html 源码文件，采用 HTML 编程语言编写，共 418 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?