zhptrs.f.html

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

HTML
426 行
字号
12
</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="ZLACGV.203"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL ZGEMV( <span class="string">'Conjugate transpose'</span>, K-1, NRHS, -ONE, B,
     $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
               CALL <a name="ZLACGV.206"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</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 ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
            KC = KC + K
            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="ZLACGV.224"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL ZGEMV( <span class="string">'Conjugate transpose'</span>, K-1, NRHS, -ONE, B,
     $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
               CALL <a name="ZLACGV.227"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLACGV.229"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K+1, 1 ), LDB )
               CALL ZGEMV( <span class="string">'Conjugate transpose'</span>, K-1, NRHS, -ONE, B,
     $                     LDB, AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
               CALL <a name="ZLACGV.232"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</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 ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
            KC = KC + 2*K + 1
            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
         KC = 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 ZSWAP( 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 ZGERU( N-K, NRHS, -ONE, AP( KC+1 ), 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 = DBLE( ONE ) / DBLE( AP( KC ) )
            CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
            KC = KC + N - K + 1
            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 ZSWAP( 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 ZGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
     $                     LDB, B( K+2, 1 ), LDB )
               CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 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 = AP( KC+1 )
            AKM1 = AP( KC ) / DCONJG( AKM1K )
            AK = AP( KC+N-K+1 ) / AKM1K
            DENOM = AKM1*AK - ONE
            DO 70 J = 1, NRHS
               BKM1 = B( K, J ) / DCONJG( AKM1K )
               BK = B( K+1, J ) / AKM1K
               B( K, J ) = ( AK*BKM1-BK ) / DENOM
               B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
   70       CONTINUE
            KC = KC + 2*( N-K ) + 1
            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
         KC = N*( N+1 ) / 2 + 1
   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>         KC = KC - ( N-K+1 )
         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="ZLACGV.350"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL ZGEMV( <span class="string">'Conjugate transpose'</span>, N-K, NRHS, -ONE,
     $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
     $                     B( K, 1 ), LDB )
               CALL <a name="ZLACGV.354"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</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 ZSWAP( 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="ZLACGV.371"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K, 1 ), LDB )
               CALL ZGEMV( <span class="string">'Conjugate transpose'</span>, N-K, NRHS, -ONE,
     $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
     $                     B( K, 1 ), LDB )
               CALL <a name="ZLACGV.375"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLACGV.377"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NRHS, B( K-1, 1 ), LDB )
               CALL ZGEMV( <span class="string">'Conjugate transpose'</span>, N-K, NRHS, -ONE,
     $                     B( K+1, 1 ), LDB, AP( KC-( N-K ) ), 1, ONE,
     $                     B( K-1, 1 ), LDB )
               CALL <a name="ZLACGV.381"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</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 ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
            KC = KC - ( N-K+2 )
            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="ZHPTRS.399"></a><a href="zhptrs.f.html#ZHPTRS.1">ZHPTRS</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
12

zhptrs.f.html - 源码说明

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

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

⌨️ 快捷键说明

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