ztgsy2.f.html

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

HTML
386 行
字号
12
      INTEGER            I, IERR, J, K
      DOUBLE PRECISION   SCALOC
      COMPLEX*16         ALPHA
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            IPIV( LDZ ), JPIV( LDZ )
      COMPLEX*16         RHS( LDZ ), Z( LDZ, LDZ )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.184"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      EXTERNAL           <a name="LSAME.185"></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="XERBLA.188"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, ZAXPY, <a name="ZGESC2.188"></a><a href="zgesc2.f.html#ZGESC2.1">ZGESC2</a>, <a name="ZGETC2.188"></a><a href="zgetc2.f.html#ZGETC2.1">ZGETC2</a>, <a name="ZLATDF.188"></a><a href="zlatdf.f.html#ZLATDF.1">ZLATDF</a>, ZSCAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          DCMPLX, DCONJG, MAX
<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">     Decode and test input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IERR = 0
      NOTRAN = <a name="LSAME.199"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> )
      IF( .NOT.NOTRAN .AND. .NOT.<a name="LSAME.200"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'C'</span> ) ) THEN
         INFO = -1
      ELSE IF( NOTRAN ) THEN
         IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
            INFO = -2
         END IF
      END IF
      IF( INFO.EQ.0 ) THEN
         IF( M.LE.0 ) THEN
            INFO = -3
         ELSE IF( N.LE.0 ) THEN
            INFO = -4
         ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
            INFO = -5
         ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
            INFO = -8
         ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
            INFO = -10
         ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
            INFO = -12
         ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
            INFO = -14
         ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
            INFO = -16
         END IF
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.227"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZTGSY2.227"></a><a href="ztgsy2.f.html#ZTGSY2.1">ZTGSY2</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( NOTRAN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve (I, J) - system
</span><span class="comment">*</span><span class="comment">           A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
</span><span class="comment">*</span><span class="comment">           D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
</span><span class="comment">*</span><span class="comment">        for I = M, M - 1, ..., 1; J = 1, 2, ..., N
</span><span class="comment">*</span><span class="comment">
</span>         SCALE = ONE
         SCALOC = ONE
         DO 30 J = 1, N
            DO 20 I = M, 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Build 2 by 2 system
</span><span class="comment">*</span><span class="comment">
</span>               Z( 1, 1 ) = A( I, I )
               Z( 2, 1 ) = D( I, I )
               Z( 1, 2 ) = -B( J, J )
               Z( 2, 2 ) = -E( J, J )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Set up right hand side(s)
</span><span class="comment">*</span><span class="comment">
</span>               RHS( 1 ) = C( I, J )
               RHS( 2 ) = F( I, J )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve Z * x = RHS
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZGETC2.257"></a><a href="zgetc2.f.html#ZGETC2.1">ZGETC2</a>( LDZ, Z, LDZ, IPIV, JPIV, IERR )
               IF( IERR.GT.0 )
     $            INFO = IERR
               IF( IJOB.EQ.0 ) THEN
                  CALL <a name="ZGESC2.261"></a><a href="zgesc2.f.html#ZGESC2.1">ZGESC2</a>( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
                  IF( SCALOC.NE.ONE ) THEN
                     DO 10 K = 1, N
                        CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
     $                              C( 1, K ), 1 )
                        CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
     $                              F( 1, K ), 1 )
   10                CONTINUE
                     SCALE = SCALE*SCALOC
                  END IF
               ELSE
                  CALL <a name="ZLATDF.272"></a><a href="zlatdf.f.html#ZLATDF.1">ZLATDF</a>( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL,
     $                         IPIV, JPIV )
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Unpack solution vector(s)
</span><span class="comment">*</span><span class="comment">
</span>               C( I, J ) = RHS( 1 )
               F( I, J ) = RHS( 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Substitute R(I, J) and L(I, J) into remaining equation.
</span><span class="comment">*</span><span class="comment">
</span>               IF( I.GT.1 ) THEN
                  ALPHA = -RHS( 1 )
                  CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 )
                  CALL ZAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 )
               END IF
               IF( J.LT.N ) THEN
                  CALL ZAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB,
     $                        C( I, J+1 ), LDC )
                  CALL ZAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE,
     $                        F( I, J+1 ), LDF )
               END IF
<span class="comment">*</span><span class="comment">
</span>   20       CONTINUE
   30    CONTINUE
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve transposed (I, J) - system:
</span><span class="comment">*</span><span class="comment">           A(I, I)' * R(I, J) + D(I, I)' * L(J, J) = C(I, J)
</span><span class="comment">*</span><span class="comment">           R(I, I) * B(J, J) + L(I, J) * E(J, J)   = -F(I, J)
</span><span class="comment">*</span><span class="comment">        for I = 1, 2, ..., M, J = N, N - 1, ..., 1
</span><span class="comment">*</span><span class="comment">
</span>         SCALE = ONE
         SCALOC = ONE
         DO 80 I = 1, M
            DO 70 J = N, 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Build 2 by 2 system Z'
</span><span class="comment">*</span><span class="comment">
</span>               Z( 1, 1 ) = DCONJG( A( I, I ) )
               Z( 2, 1 ) = -DCONJG( B( J, J ) )
               Z( 1, 2 ) = DCONJG( D( I, I ) )
               Z( 2, 2 ) = -DCONJG( E( J, J ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Set up right hand side(s)
</span><span class="comment">*</span><span class="comment">
</span>               RHS( 1 ) = C( I, J )
               RHS( 2 ) = F( I, J )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Solve Z' * x = RHS
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZGETC2.324"></a><a href="zgetc2.f.html#ZGETC2.1">ZGETC2</a>( LDZ, Z, LDZ, IPIV, JPIV, IERR )
               IF( IERR.GT.0 )
     $            INFO = IERR
               CALL <a name="ZGESC2.327"></a><a href="zgesc2.f.html#ZGESC2.1">ZGESC2</a>( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
               IF( SCALOC.NE.ONE ) THEN
                  DO 40 K = 1, N
                     CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
     $                           1 )
                     CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
     $                           1 )
   40             CONTINUE
                  SCALE = SCALE*SCALOC
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Unpack solution vector(s)
</span><span class="comment">*</span><span class="comment">
</span>               C( I, J ) = RHS( 1 )
               F( I, J ) = RHS( 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Substitute R(I, J) and L(I, J) into remaining equation.
</span><span class="comment">*</span><span class="comment">
</span>               DO 50 K = 1, J - 1
                  F( I, K ) = F( I, K ) + RHS( 1 )*DCONJG( B( K, J ) ) +
     $                        RHS( 2 )*DCONJG( E( K, J ) )
   50          CONTINUE
               DO 60 K = I + 1, M
                  C( K, J ) = C( K, J ) - DCONJG( A( I, K ) )*RHS( 1 ) -
     $                        DCONJG( D( I, K ) )*RHS( 2 )
   60          CONTINUE
<span class="comment">*</span><span class="comment">
</span>   70       CONTINUE
   80    CONTINUE
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZTGSY2.359"></a><a href="ztgsy2.f.html#ZTGSY2.1">ZTGSY2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
12

ztgsy2.f.html - 源码说明

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

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

⌨️ 快捷键说明

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