LAPACK: SRC/ctrsyl.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
00002      $                   LDC, SCALE, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          TRANA, TRANB
00011       INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
00012       REAL               SCALE
00013 *     ..
00014 *     .. Array Arguments ..
00015       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CTRSYL solves the complex Sylvester matrix equation:
00022 *
00023 *     op(A)*X + X*op(B) = scale*C or
00024 *     op(A)*X - X*op(B) = scale*C,
00025 *
00026 *  where op(A) = A or A**H, and A and B are both upper triangular. A is
00027 *  M-by-M and B is N-by-N; the right hand side C and the solution X are
00028 *  M-by-N; and scale is an output scale factor, set <= 1 to avoid
00029 *  overflow in X.
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  TRANA   (input) CHARACTER*1
00035 *          Specifies the option op(A):
00036 *          = 'N': op(A) = A    (No transpose)
00037 *          = 'C': op(A) = A**H (Conjugate transpose)
00038 *
00039 *  TRANB   (input) CHARACTER*1
00040 *          Specifies the option op(B):
00041 *          = 'N': op(B) = B    (No transpose)
00042 *          = 'C': op(B) = B**H (Conjugate transpose)
00043 *
00044 *  ISGN    (input) INTEGER
00045 *          Specifies the sign in the equation:
00046 *          = +1: solve op(A)*X + X*op(B) = scale*C
00047 *          = -1: solve op(A)*X - X*op(B) = scale*C
00048 *
00049 *  M       (input) INTEGER
00050 *          The order of the matrix A, and the number of rows in the
00051 *          matrices X and C. M >= 0.
00052 *
00053 *  N       (input) INTEGER
00054 *          The order of the matrix B, and the number of columns in the
00055 *          matrices X and C. N >= 0.
00056 *
00057 *  A       (input) COMPLEX array, dimension (LDA,M)
00058 *          The upper triangular matrix A.
00059 *
00060 *  LDA     (input) INTEGER
00061 *          The leading dimension of the array A. LDA >= max(1,M).
00062 *
00063 *  B       (input) COMPLEX array, dimension (LDB,N)
00064 *          The upper triangular matrix B.
00065 *
00066 *  LDB     (input) INTEGER
00067 *          The leading dimension of the array B. LDB >= max(1,N).
00068 *
00069 *  C       (input/output) COMPLEX array, dimension (LDC,N)
00070 *          On entry, the M-by-N right hand side matrix C.
00071 *          On exit, C is overwritten by the solution matrix X.
00072 *
00073 *  LDC     (input) INTEGER
00074 *          The leading dimension of the array C. LDC >= max(1,M)
00075 *
00076 *  SCALE   (output) REAL
00077 *          The scale factor, scale, set <= 1 to avoid overflow in X.
00078 *
00079 *  INFO    (output) INTEGER
00080 *          = 0: successful exit
00081 *          < 0: if INFO = -i, the i-th argument had an illegal value
00082 *          = 1: A and B have common or very close eigenvalues; perturbed
00083 *               values were used to solve the equation (but the matrices
00084 *               A and B are unchanged).
00085 *
00086 *  =====================================================================
00087 *
00088 *     .. Parameters ..
00089       REAL               ONE
00090       PARAMETER          ( ONE = 1.0E+0 )
00091 *     ..
00092 *     .. Local Scalars ..
00093       LOGICAL            NOTRNA, NOTRNB
00094       INTEGER            J, K, L
00095       REAL               BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
00096      $                   SMLNUM
00097       COMPLEX            A11, SUML, SUMR, VEC, X11
00098 *     ..
00099 *     .. Local Arrays ..
00100       REAL               DUM( 1 )
00101 *     ..
00102 *     .. External Functions ..
00103       LOGICAL            LSAME
00104       REAL               CLANGE, SLAMCH
00105       COMPLEX            CDOTC, CDOTU, CLADIV
00106       EXTERNAL           LSAME, CLANGE, SLAMCH, CDOTC, CDOTU, CLADIV
00107 *     ..
00108 *     .. External Subroutines ..
00109       EXTERNAL           CSSCAL, SLABAD, XERBLA
00110 *     ..
00111 *     .. Intrinsic Functions ..
00112       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
00113 *     ..
00114 *     .. Executable Statements ..
00115 *
00116 *     Decode and Test input parameters
00117 *
00118       NOTRNA = LSAME( TRANA, 'N' )
00119       NOTRNB = LSAME( TRANB, 'N' )
00120 *
00121       INFO = 0
00122       IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'C' ) ) THEN
00123          INFO = -1
00124       ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'C' ) ) THEN
00125          INFO = -2
00126       ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
00127          INFO = -3
00128       ELSE IF( M.LT.0 ) THEN
00129          INFO = -4
00130       ELSE IF( N.LT.0 ) THEN
00131          INFO = -5
00132       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00133          INFO = -7
00134       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00135          INFO = -9
00136       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
00137          INFO = -11
00138       END IF
00139       IF( INFO.NE.0 ) THEN
00140          CALL XERBLA( 'CTRSYL', -INFO )
00141          RETURN
00142       END IF
00143 *
00144 *     Quick return if possible
00145 *
00146       SCALE = ONE
00147       IF( M.EQ.0 .OR. N.EQ.0 )
00148      $   RETURN
00149 *
00150 *     Set constants to control overflow
00151 *
00152       EPS = SLAMCH( 'P' )
00153       SMLNUM = SLAMCH( 'S' )
00154       BIGNUM = ONE / SMLNUM
00155       CALL SLABAD( SMLNUM, BIGNUM )
00156       SMLNUM = SMLNUM*REAL( M*N ) / EPS
00157       BIGNUM = ONE / SMLNUM
00158       SMIN = MAX( SMLNUM, EPS*CLANGE( 'M', M, M, A, LDA, DUM ),
00159      $       EPS*CLANGE( 'M', N, N, B, LDB, DUM ) )
00160       SGN = ISGN
00161 *
00162       IF( NOTRNA .AND. NOTRNB ) THEN
00163 *
00164 *        Solve    A*X + ISGN*X*B = scale*C.
00165 *
00166 *        The (K,L)th block of X is determined starting from
00167 *        bottom-left corner column by column by
00168 *
00169 *            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
00170 *
00171 *        Where
00172 *                    M                        L-1
00173 *          R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)].
00174 *                  I=K+1                      J=1
00175 *
00176          DO 30 L = 1, N
00177             DO 20 K = M, 1, -1
00178 *
00179                SUML = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
00180      $                C( MIN( K+1, M ), L ), 1 )
00181                SUMR = CDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
00182                VEC = C( K, L ) - ( SUML+SGN*SUMR )
00183 *
00184                SCALOC = ONE
00185                A11 = A( K, K ) + SGN*B( L, L )
00186                DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
00187                IF( DA11.LE.SMIN ) THEN
00188                   A11 = SMIN
00189                   DA11 = SMIN
00190                   INFO = 1
00191                END IF
00192                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
00193                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
00194                   IF( DB.GT.BIGNUM*DA11 )
00195      $               SCALOC = ONE / DB
00196                END IF
00197                X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
00198 *
00199                IF( SCALOC.NE.ONE ) THEN
00200                   DO 10 J = 1, N
00201                      CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
00202    10             CONTINUE
00203                   SCALE = SCALE*SCALOC
00204                END IF
00205                C( K, L ) = X11
00206 *
00207    20       CONTINUE
00208    30    CONTINUE
00209 *
00210       ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
00211 *
00212 *        Solve    A' *X + ISGN*X*B = scale*C.
00213 *
00214 *        The (K,L)th block of X is determined starting from
00215 *        upper-left corner column by column by
00216 *
00217 *            A'(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
00218 *
00219 *        Where
00220 *                   K-1                         L-1
00221 *          R(K,L) = SUM [A'(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]
00222 *                   I=1                         J=1
00223 *
00224          DO 60 L = 1, N
00225             DO 50 K = 1, M
00226 *
00227                SUML = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
00228                SUMR = CDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
00229                VEC = C( K, L ) - ( SUML+SGN*SUMR )
00230 *
00231                SCALOC = ONE
00232                A11 = CONJG( A( K, K ) ) + SGN*B( L, L )
00233                DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
00234                IF( DA11.LE.SMIN ) THEN
00235                   A11 = SMIN
00236                   DA11 = SMIN
00237                   INFO = 1
00238                END IF
00239                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
00240                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
00241                   IF( DB.GT.BIGNUM*DA11 )
00242      $               SCALOC = ONE / DB
00243                END IF
00244 *
00245                X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
00246 *
00247                IF( SCALOC.NE.ONE ) THEN
00248                   DO 40 J = 1, N
00249                      CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
00250    40             CONTINUE
00251                   SCALE = SCALE*SCALOC
00252                END IF
00253                C( K, L ) = X11
00254 *
00255    50       CONTINUE
00256    60    CONTINUE
00257 *
00258       ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
00259 *
00260 *        Solve    A'*X + ISGN*X*B' = C.
00261 *
00262 *        The (K,L)th block of X is determined starting from
00263 *        upper-right corner column by column by
00264 *
00265 *            A'(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
00266 *
00267 *        Where
00268 *                    K-1
00269 *           R(K,L) = SUM [A'(I,K)*X(I,L)] +
00270 *                    I=1
00271 *                           N
00272 *                     ISGN*SUM [X(K,J)*B'(L,J)].
00273 *                          J=L+1
00274 *
00275          DO 90 L = N, 1, -1
00276             DO 80 K = 1, M
00277 *
00278                SUML = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
00279                SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
00280      $                B( L, MIN( L+1, N ) ), LDB )
00281                VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) )
00282 *
00283                SCALOC = ONE
00284                A11 = CONJG( A( K, K )+SGN*B( L, L ) )
00285                DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
00286                IF( DA11.LE.SMIN ) THEN
00287                   A11 = SMIN
00288                   DA11 = SMIN
00289                   INFO = 1
00290                END IF
00291                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
00292                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
00293                   IF( DB.GT.BIGNUM*DA11 )
00294      $               SCALOC = ONE / DB
00295                END IF
00296 *
00297                X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
00298 *
00299                IF( SCALOC.NE.ONE ) THEN
00300                   DO 70 J = 1, N
00301                      CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
00302    70             CONTINUE
00303                   SCALE = SCALE*SCALOC
00304                END IF
00305                C( K, L ) = X11
00306 *
00307    80       CONTINUE
00308    90    CONTINUE
00309 *
00310       ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
00311 *
00312 *        Solve    A*X + ISGN*X*B' = C.
00313 *
00314 *        The (K,L)th block of X is determined starting from
00315 *        bottom-left corner column by column by
00316 *
00317 *           A(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
00318 *
00319 *        Where
00320 *                    M                          N
00321 *          R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B'(L,J)]
00322 *                  I=K+1                      J=L+1
00323 *
00324          DO 120 L = N, 1, -1
00325             DO 110 K = M, 1, -1
00326 *
00327                SUML = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
00328      $                C( MIN( K+1, M ), L ), 1 )
00329                SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
00330      $                B( L, MIN( L+1, N ) ), LDB )
00331                VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) )
00332 *
00333                SCALOC = ONE
00334                A11 = A( K, K ) + SGN*CONJG( B( L, L ) )
00335                DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
00336                IF( DA11.LE.SMIN ) THEN
00337                   A11 = SMIN
00338                   DA11 = SMIN
00339                   INFO = 1
00340                END IF
00341                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) )
00342                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
00343                   IF( DB.GT.BIGNUM*DA11 )
00344      $               SCALOC = ONE / DB
00345                END IF
00346 *
00347                X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
00348 *
00349                IF( SCALOC.NE.ONE ) THEN
00350                   DO 100 J = 1, N
00351                      CALL CSSCAL( M, SCALOC, C( 1, J ), 1 )
00352   100             CONTINUE
00353                   SCALE = SCALE*SCALOC
00354                END IF
00355                C( K, L ) = X11
00356 *
00357   110       CONTINUE
00358   120    CONTINUE
00359 *
00360       END IF
00361 *
00362       RETURN
00363 *
00364 *     End of CTRSYL
00365 *
00366       END