ctrsyl.f.html

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

<span class="comment">*</span><span class="comment">
</span>               SUML = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
     $                C( MIN( K+1, M ), L ), 1 )
               SUMR = CDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
               VEC = C( K, L ) - ( SUML+SGN*SUMR )
<span class="comment">*</span><span class="comment">
</span>               SCALOC = ONE
               A11 = A( K, K ) + SGN*B( L, L )
               DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
               IF( DA11.LE.SMIN ) THEN
                  A11 = SMIN
                  DA11 = SMIN
                  INFO = 1
               END IF
               DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
               IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
                  IF( DB.GT.BIGNUM*DA11 )
     $               SCALOC = ONE / DB
               END IF
               X11 = <a name="CLADIV.196"></a><a href="cladiv.f.html#CLADIV.1">CLADIV</a>( VEC*CMPLX( SCALOC ), A11 )
<span class="comment">*</span><span class="comment">
</span>               IF( SCALOC.NE.ONE ) THEN
                  DO 10 J = 1, N
                     CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
   10             CONTINUE
                  SCALE = SCALE*SCALOC
               END IF
               C( K, L ) = X11
<span class="comment">*</span><span class="comment">
</span>   20       CONTINUE
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve    A' *X + ISGN*X*B = scale*C.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        The (K,L)th block of X is determined starting from
</span><span class="comment">*</span><span class="comment">        upper-left corner column by column by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            A'(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Where
</span><span class="comment">*</span><span class="comment">                   K-1                         L-1
</span><span class="comment">*</span><span class="comment">          R(K,L) = SUM [A'(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]
</span><span class="comment">*</span><span class="comment">                   I=1                         J=1
</span><span class="comment">*</span><span class="comment">
</span>         DO 60 L = 1, N
            DO 50 K = 1, M
<span class="comment">*</span><span class="comment">
</span>               SUML = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
               SUMR = CDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
               VEC = C( K, L ) - ( SUML+SGN*SUMR )
<span class="comment">*</span><span class="comment">
</span>               SCALOC = ONE
               A11 = CONJG( A( K, K ) ) + SGN*B( L, L )
               DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
               IF( DA11.LE.SMIN ) THEN
                  A11 = SMIN
                  DA11 = SMIN
                  INFO = 1
               END IF
               DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
               IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
                  IF( DB.GT.BIGNUM*DA11 )
     $               SCALOC = ONE / DB
               END IF
<span class="comment">*</span><span class="comment">
</span>               X11 = <a name="CLADIV.244"></a><a href="cladiv.f.html#CLADIV.1">CLADIV</a>( VEC*CMPLX( SCALOC ), A11 )
<span class="comment">*</span><span class="comment">
</span>               IF( SCALOC.NE.ONE ) THEN
                  DO 40 J = 1, N
                     CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
   40             CONTINUE
                  SCALE = SCALE*SCALOC
               END IF
               C( K, L ) = X11
<span class="comment">*</span><span class="comment">
</span>   50       CONTINUE
   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve    A'*X + ISGN*X*B' = C.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        The (K,L)th block of X is determined starting from
</span><span class="comment">*</span><span class="comment">        upper-right corner column by column by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            A'(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Where
</span><span class="comment">*</span><span class="comment">                    K-1
</span><span class="comment">*</span><span class="comment">           R(K,L) = SUM [A'(I,K)*X(I,L)] +
</span><span class="comment">*</span><span class="comment">                    I=1
</span><span class="comment">*</span><span class="comment">                           N
</span><span class="comment">*</span><span class="comment">                     ISGN*SUM [X(K,J)*B'(L,J)].
</span><span class="comment">*</span><span class="comment">                          J=L+1
</span><span class="comment">*</span><span class="comment">
</span>         DO 90 L = N, 1, -1
            DO 80 K = 1, M
<span class="comment">*</span><span class="comment">
</span>               SUML = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
               SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
     $                B( L, MIN( L+1, N ) ), LDB )
               VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) )
<span class="comment">*</span><span class="comment">
</span>               SCALOC = ONE
               A11 = CONJG( A( K, K )+SGN*B( L, L ) )
               DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
               IF( DA11.LE.SMIN ) THEN
                  A11 = SMIN
                  DA11 = SMIN
                  INFO = 1
               END IF
               DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
               IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
                  IF( DB.GT.BIGNUM*DA11 )
     $               SCALOC = ONE / DB
               END IF
<span class="comment">*</span><span class="comment">
</span>               X11 = <a name="CLADIV.296"></a><a href="cladiv.f.html#CLADIV.1">CLADIV</a>( VEC*CMPLX( SCALOC ), A11 )
<span class="comment">*</span><span class="comment">
</span>               IF( SCALOC.NE.ONE ) THEN
                  DO 70 J = 1, N
                     CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
   70             CONTINUE
                  SCALE = SCALE*SCALOC
               END IF
               C( K, L ) = X11
<span class="comment">*</span><span class="comment">
</span>   80       CONTINUE
   90    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve    A*X + ISGN*X*B' = C.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        The (K,L)th block of X is determined starting from
</span><span class="comment">*</span><span class="comment">        bottom-left corner column by column by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           A(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Where
</span><span class="comment">*</span><span class="comment">                    M                          N
</span><span class="comment">*</span><span class="comment">          R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B'(L,J)]
</span><span class="comment">*</span><span class="comment">                  I=K+1                      J=L+1
</span><span class="comment">*</span><span class="comment">
</span>         DO 120 L = N, 1, -1
            DO 110 K = M, 1, -1
<span class="comment">*</span><span class="comment">
</span>               SUML = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
     $                C( MIN( K+1, M ), L ), 1 )
               SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
     $                B( L, MIN( L+1, N ) ), LDB )
               VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) )
<span class="comment">*</span><span class="comment">
</span>               SCALOC = ONE
               A11 = A( K, K ) + SGN*CONJG( B( L, L ) )
               DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
               IF( DA11.LE.SMIN ) THEN
                  A11 = SMIN
                  DA11 = SMIN
                  INFO = 1
               END IF
               DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
               IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
                  IF( DB.GT.BIGNUM*DA11 )
     $               SCALOC = ONE / DB
               END IF
<span class="comment">*</span><span class="comment">
</span>               X11 = <a name="CLADIV.346"></a><a href="cladiv.f.html#CLADIV.1">CLADIV</a>( VEC*CMPLX( SCALOC ), A11 )
<span class="comment">*</span><span class="comment">
</span>               IF( SCALOC.NE.ONE ) THEN
                  DO 100 J = 1, N
                     CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
  100             CONTINUE
                  SCALE = SCALE*SCALOC
               END IF
               C( K, L ) = X11
<span class="comment">*</span><span class="comment">
</span>  110       CONTINUE
  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      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="CTRSYL.363"></a><a href="ctrsyl.f.html#CTRSYL.1">CTRSYL</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>