|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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
1.7.3 
