1 SUBROUTINE CTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 2 $ LDC, SCALE, INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER TRANA, TRANB 11 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 12 REAL SCALE 13* .. 14* .. Array Arguments .. 15 COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* CTRSYL solves the complex Sylvester matrix equation: 22* 23* op(A)*X + X*op(B) = scale*C or 24* op(A)*X - X*op(B) = scale*C, 25* 26* where op(A) = A or A**H, and A and B are both upper triangular. A is 27* M-by-M and B is N-by-N; the right hand side C and the solution X are 28* M-by-N; and scale is an output scale factor, set <= 1 to avoid 29* overflow in X. 30* 31* Arguments 32* ========= 33* 34* TRANA (input) CHARACTER*1 35* Specifies the option op(A): 36* = 'N': op(A) = A (No transpose) 37* = 'C': op(A) = A**H (Conjugate transpose) 38* 39* TRANB (input) CHARACTER*1 40* Specifies the option op(B): 41* = 'N': op(B) = B (No transpose) 42* = 'C': op(B) = B**H (Conjugate transpose) 43* 44* ISGN (input) INTEGER 45* Specifies the sign in the equation: 46* = +1: solve op(A)*X + X*op(B) = scale*C 47* = -1: solve op(A)*X - X*op(B) = scale*C 48* 49* M (input) INTEGER 50* The order of the matrix A, and the number of rows in the 51* matrices X and C. M >= 0. 52* 53* N (input) INTEGER 54* The order of the matrix B, and the number of columns in the 55* matrices X and C. N >= 0. 56* 57* A (input) COMPLEX array, dimension (LDA,M) 58* The upper triangular matrix A. 59* 60* LDA (input) INTEGER 61* The leading dimension of the array A. LDA >= max(1,M). 62* 63* B (input) COMPLEX array, dimension (LDB,N) 64* The upper triangular matrix B. 65* 66* LDB (input) INTEGER 67* The leading dimension of the array B. LDB >= max(1,N). 68* 69* C (input/output) COMPLEX array, dimension (LDC,N) 70* On entry, the M-by-N right hand side matrix C. 71* On exit, C is overwritten by the solution matrix X. 72* 73* LDC (input) INTEGER 74* The leading dimension of the array C. LDC >= max(1,M) 75* 76* SCALE (output) REAL 77* The scale factor, scale, set <= 1 to avoid overflow in X. 78* 79* INFO (output) INTEGER 80* = 0: successful exit 81* < 0: if INFO = -i, the i-th argument had an illegal value 82* = 1: A and B have common or very close eigenvalues; perturbed 83* values were used to solve the equation (but the matrices 84* A and B are unchanged). 85* 86* ===================================================================== 87* 88* .. Parameters .. 89 REAL ONE 90 PARAMETER ( ONE = 1.0E+0 ) 91* .. 92* .. Local Scalars .. 93 LOGICAL NOTRNA, NOTRNB 94 INTEGER J, K, L 95 REAL BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 96 $ SMLNUM 97 COMPLEX A11, SUML, SUMR, VEC, X11 98* .. 99* .. Local Arrays .. 100 REAL DUM( 1 ) 101* .. 102* .. External Functions .. 103 LOGICAL LSAME 104 REAL CLANGE, SLAMCH 105 COMPLEX CDOTC, CDOTU, CLADIV 106 EXTERNAL LSAME, CLANGE, SLAMCH, CDOTC, CDOTU, CLADIV 107* .. 108* .. External Subroutines .. 109 EXTERNAL CSSCAL, SLABAD, XERBLA 110* .. 111* .. Intrinsic Functions .. 112 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL 113* .. 114* .. Executable Statements .. 115* 116* Decode and Test input parameters 117* 118 NOTRNA = LSAME( TRANA, 'N' ) 119 NOTRNB = LSAME( TRANB, 'N' ) 120* 121 INFO = 0 122 IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT. 123 $ LSAME( TRANA, 'C' ) ) THEN 124 INFO = -1 125 ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT. 126 $ LSAME( TRANB, 'C' ) ) THEN 127 INFO = -2 128 ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN 129 INFO = -3 130 ELSE IF( M.LT.0 ) THEN 131 INFO = -4 132 ELSE IF( N.LT.0 ) THEN 133 INFO = -5 134 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 135 INFO = -7 136 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 137 INFO = -9 138 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 139 INFO = -11 140 END IF 141 IF( INFO.NE.0 ) THEN 142 CALL XERBLA( 'CTRSYL', -INFO ) 143 RETURN 144 END IF 145* 146* Quick return if possible 147* 148 IF( M.EQ.0 .OR. N.EQ.0 ) 149 $ RETURN 150* 151* Set constants to control overflow 152* 153 EPS = SLAMCH( 'P' ) 154 SMLNUM = SLAMCH( 'S' ) 155 BIGNUM = ONE / SMLNUM 156 CALL SLABAD( SMLNUM, BIGNUM ) 157 SMLNUM = SMLNUM*REAL( M*N ) / EPS 158 BIGNUM = ONE / SMLNUM 159 SMIN = MAX( SMLNUM, EPS*CLANGE( 'M', M, M, A, LDA, DUM ), 160 $ EPS*CLANGE( 'M', N, N, B, LDB, DUM ) ) 161 SCALE = ONE 162 SGN = ISGN 163* 164 IF( NOTRNA .AND. NOTRNB ) THEN 165* 166* Solve A*X + ISGN*X*B = scale*C. 167* 168* The (K,L)th block of X is determined starting from 169* bottom-left corner column by column by 170* 171* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 172* 173* Where 174* M L-1 175* R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]. 176* I=K+1 J=1 177* 178 DO 30 L = 1, N 179 DO 20 K = M, 1, -1 180* 181 SUML = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA, 182 $ C( MIN( K+1, M ), L ), 1 ) 183 SUMR = CDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 ) 184 VEC = C( K, L ) - ( SUML+SGN*SUMR ) 185* 186 SCALOC = ONE 187 A11 = A( K, K ) + SGN*B( L, L ) 188 DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) ) 189 IF( DA11.LE.SMIN ) THEN 190 A11 = SMIN 191 DA11 = SMIN 192 INFO = 1 193 END IF 194 DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) ) 195 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 196 IF( DB.GT.BIGNUM*DA11 ) 197 $ SCALOC = ONE / DB 198 END IF 199 X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 ) 200* 201 IF( SCALOC.NE.ONE ) THEN 202 DO 10 J = 1, N 203 CALL CSSCAL( M, SCALOC, C( 1, J ), 1 ) 204 10 CONTINUE 205 SCALE = SCALE*SCALOC 206 END IF 207 C( K, L ) = X11 208* 209 20 CONTINUE 210 30 CONTINUE 211* 212 ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN 213* 214* Solve A' *X + ISGN*X*B = scale*C. 215* 216* The (K,L)th block of X is determined starting from 217* upper-left corner column by column by 218* 219* A'(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 220* 221* Where 222* K-1 L-1 223* R(K,L) = SUM [A'(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)] 224* I=1 J=1 225* 226 DO 60 L = 1, N 227 DO 50 K = 1, M 228* 229 SUML = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 ) 230 SUMR = CDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 ) 231 VEC = C( K, L ) - ( SUML+SGN*SUMR ) 232* 233 SCALOC = ONE 234 A11 = CONJG( A( K, K ) ) + SGN*B( L, L ) 235 DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) ) 236 IF( DA11.LE.SMIN ) THEN 237 A11 = SMIN 238 DA11 = SMIN 239 INFO = 1 240 END IF 241 DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) ) 242 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 243 IF( DB.GT.BIGNUM*DA11 ) 244 $ SCALOC = ONE / DB 245 END IF 246* 247 X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 ) 248* 249 IF( SCALOC.NE.ONE ) THEN 250 DO 40 J = 1, N 251 CALL CSSCAL( M, SCALOC, C( 1, J ), 1 ) 252 40 CONTINUE 253 SCALE = SCALE*SCALOC 254 END IF 255 C( K, L ) = X11 256* 257 50 CONTINUE 258 60 CONTINUE 259* 260 ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN 261* 262* Solve A'*X + ISGN*X*B' = C. 263* 264* The (K,L)th block of X is determined starting from 265* upper-right corner column by column by 266* 267* A'(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L) 268* 269* Where 270* K-1 271* R(K,L) = SUM [A'(I,K)*X(I,L)] + 272* I=1 273* N 274* ISGN*SUM [X(K,J)*B'(L,J)]. 275* J=L+1 276* 277 DO 90 L = N, 1, -1 278 DO 80 K = 1, M 279* 280 SUML = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 ) 281 SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC, 282 $ B( L, MIN( L+1, N ) ), LDB ) 283 VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) ) 284* 285 SCALOC = ONE 286 A11 = CONJG( A( K, K )+SGN*B( L, L ) ) 287 DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) ) 288 IF( DA11.LE.SMIN ) THEN 289 A11 = SMIN 290 DA11 = SMIN 291 INFO = 1 292 END IF 293 DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) ) 294 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 295 IF( DB.GT.BIGNUM*DA11 ) 296 $ SCALOC = ONE / DB 297 END IF 298* 299 X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 ) 300* 301 IF( SCALOC.NE.ONE ) THEN 302 DO 70 J = 1, N 303 CALL CSSCAL( M, SCALOC, C( 1, J ), 1 ) 304 70 CONTINUE 305 SCALE = SCALE*SCALOC 306 END IF 307 C( K, L ) = X11 308* 309 80 CONTINUE 310 90 CONTINUE 311* 312 ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN 313* 314* Solve A*X + ISGN*X*B' = C. 315* 316* The (K,L)th block of X is determined starting from 317* bottom-left corner column by column by 318* 319* A(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L) 320* 321* Where 322* M N 323* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B'(L,J)] 324* I=K+1 J=L+1 325* 326 DO 120 L = N, 1, -1 327 DO 110 K = M, 1, -1 328* 329 SUML = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA, 330 $ C( MIN( K+1, M ), L ), 1 ) 331 SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC, 332 $ B( L, MIN( L+1, N ) ), LDB ) 333 VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) ) 334* 335 SCALOC = ONE 336 A11 = A( K, K ) + SGN*CONJG( B( L, L ) ) 337 DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) ) 338 IF( DA11.LE.SMIN ) THEN 339 A11 = SMIN 340 DA11 = SMIN 341 INFO = 1 342 END IF 343 DB = ABS( REAL( VEC ) ) + ABS( AIMAG( VEC ) ) 344 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 345 IF( DB.GT.BIGNUM*DA11 ) 346 $ SCALOC = ONE / DB 347 END IF 348* 349 X11 = CLADIV( VEC*CMPLX( SCALOC ), A11 ) 350* 351 IF( SCALOC.NE.ONE ) THEN 352 DO 100 J = 1, N 353 CALL CSSCAL( M, SCALOC, C( 1, J ), 1 ) 354 100 CONTINUE 355 SCALE = SCALE*SCALOC 356 END IF 357 C( K, L ) = X11 358* 359 110 CONTINUE 360 120 CONTINUE 361* 362 END IF 363* 364 RETURN 365* 366* End of CTRSYL 367* 368 END 369