![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 2: $ LDC, SCALE, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: CHARACTER TRANA, TRANB 11: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 12: DOUBLE PRECISION SCALE 13: * .. 14: * .. Array Arguments .. 15: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * ZTRSYL 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*16 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*16 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*16 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) DOUBLE PRECISION 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: DOUBLE PRECISION ONE 90: PARAMETER ( ONE = 1.0D+0 ) 91: * .. 92: * .. Local Scalars .. 93: LOGICAL NOTRNA, NOTRNB 94: INTEGER J, K, L 95: DOUBLE PRECISION BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 96: $ SMLNUM 97: COMPLEX*16 A11, SUML, SUMR, VEC, X11 98: * .. 99: * .. Local Arrays .. 100: DOUBLE PRECISION DUM( 1 ) 101: * .. 102: * .. External Functions .. 103: LOGICAL LSAME 104: DOUBLE PRECISION DLAMCH, ZLANGE 105: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV 106: EXTERNAL LSAME, DLAMCH, ZLANGE, ZDOTC, ZDOTU, ZLADIV 107: * .. 108: * .. External Subroutines .. 109: EXTERNAL DLABAD, XERBLA, ZDSCAL 110: * .. 111: * .. Intrinsic Functions .. 112: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN 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, 'C' ) ) THEN 123: INFO = -1 124: ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'C' ) ) THEN 125: INFO = -2 126: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN 127: INFO = -3 128: ELSE IF( M.LT.0 ) THEN 129: INFO = -4 130: ELSE IF( N.LT.0 ) THEN 131: INFO = -5 132: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 133: INFO = -7 134: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 135: INFO = -9 136: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 137: INFO = -11 138: END IF 139: IF( INFO.NE.0 ) THEN 140: CALL XERBLA( 'ZTRSYL', -INFO ) 141: RETURN 142: END IF 143: * 144: * Quick return if possible 145: * 146: SCALE = ONE 147: IF( M.EQ.0 .OR. N.EQ.0 ) 148: $ RETURN 149: * 150: * Set constants to control overflow 151: * 152: EPS = DLAMCH( 'P' ) 153: SMLNUM = DLAMCH( 'S' ) 154: BIGNUM = ONE / SMLNUM 155: CALL DLABAD( SMLNUM, BIGNUM ) 156: SMLNUM = SMLNUM*DBLE( M*N ) / EPS 157: BIGNUM = ONE / SMLNUM 158: SMIN = MAX( SMLNUM, EPS*ZLANGE( 'M', M, M, A, LDA, DUM ), 159: $ EPS*ZLANGE( 'M', N, N, B, LDB, DUM ) ) 160: SGN = ISGN 161: * 162: IF( NOTRNA .AND. NOTRNB ) THEN 163: * 164: * Solve A*X + ISGN*X*B = scale*C. 165: * 166: * The (K,L)th block of X is determined starting from 167: * bottom-left corner column by column by 168: * 169: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 170: * 171: * Where 172: * M L-1 173: * R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]. 174: * I=K+1 J=1 175: * 176: DO 30 L = 1, N 177: DO 20 K = M, 1, -1 178: * 179: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA, 180: $ C( MIN( K+1, M ), L ), 1 ) 181: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 ) 182: VEC = C( K, L ) - ( SUML+SGN*SUMR ) 183: * 184: SCALOC = ONE 185: A11 = A( K, K ) + SGN*B( L, L ) 186: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 187: IF( DA11.LE.SMIN ) THEN 188: A11 = SMIN 189: DA11 = SMIN 190: INFO = 1 191: END IF 192: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 193: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 194: IF( DB.GT.BIGNUM*DA11 ) 195: $ SCALOC = ONE / DB 196: END IF 197: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 198: * 199: IF( SCALOC.NE.ONE ) THEN 200: DO 10 J = 1, N 201: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 202: 10 CONTINUE 203: SCALE = SCALE*SCALOC 204: END IF 205: C( K, L ) = X11 206: * 207: 20 CONTINUE 208: 30 CONTINUE 209: * 210: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN 211: * 212: * Solve A' *X + ISGN*X*B = scale*C. 213: * 214: * The (K,L)th block of X is determined starting from 215: * upper-left corner column by column by 216: * 217: * A'(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 218: * 219: * Where 220: * K-1 L-1 221: * R(K,L) = SUM [A'(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)] 222: * I=1 J=1 223: * 224: DO 60 L = 1, N 225: DO 50 K = 1, M 226: * 227: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 ) 228: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 ) 229: VEC = C( K, L ) - ( SUML+SGN*SUMR ) 230: * 231: SCALOC = ONE 232: A11 = DCONJG( A( K, K ) ) + SGN*B( L, L ) 233: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 234: IF( DA11.LE.SMIN ) THEN 235: A11 = SMIN 236: DA11 = SMIN 237: INFO = 1 238: END IF 239: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 240: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 241: IF( DB.GT.BIGNUM*DA11 ) 242: $ SCALOC = ONE / DB 243: END IF 244: * 245: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 246: * 247: IF( SCALOC.NE.ONE ) THEN 248: DO 40 J = 1, N 249: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 250: 40 CONTINUE 251: SCALE = SCALE*SCALOC 252: END IF 253: C( K, L ) = X11 254: * 255: 50 CONTINUE 256: 60 CONTINUE 257: * 258: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN 259: * 260: * Solve A'*X + ISGN*X*B' = C. 261: * 262: * The (K,L)th block of X is determined starting from 263: * upper-right corner column by column by 264: * 265: * A'(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L) 266: * 267: * Where 268: * K-1 269: * R(K,L) = SUM [A'(I,K)*X(I,L)] + 270: * I=1 271: * N 272: * ISGN*SUM [X(K,J)*B'(L,J)]. 273: * J=L+1 274: * 275: DO 90 L = N, 1, -1 276: DO 80 K = 1, M 277: * 278: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 ) 279: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC, 280: $ B( L, MIN( L+1, N ) ), LDB ) 281: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) ) 282: * 283: SCALOC = ONE 284: A11 = DCONJG( A( K, K )+SGN*B( L, L ) ) 285: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 286: IF( DA11.LE.SMIN ) THEN 287: A11 = SMIN 288: DA11 = SMIN 289: INFO = 1 290: END IF 291: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 292: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 293: IF( DB.GT.BIGNUM*DA11 ) 294: $ SCALOC = ONE / DB 295: END IF 296: * 297: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 298: * 299: IF( SCALOC.NE.ONE ) THEN 300: DO 70 J = 1, N 301: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 302: 70 CONTINUE 303: SCALE = SCALE*SCALOC 304: END IF 305: C( K, L ) = X11 306: * 307: 80 CONTINUE 308: 90 CONTINUE 309: * 310: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN 311: * 312: * Solve A*X + ISGN*X*B' = C. 313: * 314: * The (K,L)th block of X is determined starting from 315: * bottom-left corner column by column by 316: * 317: * A(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L) 318: * 319: * Where 320: * M N 321: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B'(L,J)] 322: * I=K+1 J=L+1 323: * 324: DO 120 L = N, 1, -1 325: DO 110 K = M, 1, -1 326: * 327: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA, 328: $ C( MIN( K+1, M ), L ), 1 ) 329: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC, 330: $ B( L, MIN( L+1, N ) ), LDB ) 331: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) ) 332: * 333: SCALOC = ONE 334: A11 = A( K, K ) + SGN*DCONJG( B( L, L ) ) 335: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 336: IF( DA11.LE.SMIN ) THEN 337: A11 = SMIN 338: DA11 = SMIN 339: INFO = 1 340: END IF 341: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 342: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 343: IF( DB.GT.BIGNUM*DA11 ) 344: $ SCALOC = ONE / DB 345: END IF 346: * 347: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 348: * 349: IF( SCALOC.NE.ONE ) THEN 350: DO 100 J = 1, N 351: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 352: 100 CONTINUE 353: SCALE = SCALE*SCALOC 354: END IF 355: C( K, L ) = X11 356: * 357: 110 CONTINUE 358: 120 CONTINUE 359: * 360: END IF 361: * 362: RETURN 363: * 364: * End of ZTRSYL 365: * 366: END