![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTRSYL 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTRSYL + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsyl.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsyl.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsyl.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 00022 * LDC, SCALE, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER TRANA, TRANB 00026 * INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 00027 * DOUBLE PRECISION SCALE 00028 * .. 00029 * .. Array Arguments .. 00030 * COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> ZTRSYL solves the complex Sylvester matrix equation: 00040 *> 00041 *> op(A)*X + X*op(B) = scale*C or 00042 *> op(A)*X - X*op(B) = scale*C, 00043 *> 00044 *> where op(A) = A or A**H, and A and B are both upper triangular. A is 00045 *> M-by-M and B is N-by-N; the right hand side C and the solution X are 00046 *> M-by-N; and scale is an output scale factor, set <= 1 to avoid 00047 *> overflow in X. 00048 *> \endverbatim 00049 * 00050 * Arguments: 00051 * ========== 00052 * 00053 *> \param[in] TRANA 00054 *> \verbatim 00055 *> TRANA is CHARACTER*1 00056 *> Specifies the option op(A): 00057 *> = 'N': op(A) = A (No transpose) 00058 *> = 'C': op(A) = A**H (Conjugate transpose) 00059 *> \endverbatim 00060 *> 00061 *> \param[in] TRANB 00062 *> \verbatim 00063 *> TRANB is CHARACTER*1 00064 *> Specifies the option op(B): 00065 *> = 'N': op(B) = B (No transpose) 00066 *> = 'C': op(B) = B**H (Conjugate transpose) 00067 *> \endverbatim 00068 *> 00069 *> \param[in] ISGN 00070 *> \verbatim 00071 *> ISGN is INTEGER 00072 *> Specifies the sign in the equation: 00073 *> = +1: solve op(A)*X + X*op(B) = scale*C 00074 *> = -1: solve op(A)*X - X*op(B) = scale*C 00075 *> \endverbatim 00076 *> 00077 *> \param[in] M 00078 *> \verbatim 00079 *> M is INTEGER 00080 *> The order of the matrix A, and the number of rows in the 00081 *> matrices X and C. M >= 0. 00082 *> \endverbatim 00083 *> 00084 *> \param[in] N 00085 *> \verbatim 00086 *> N is INTEGER 00087 *> The order of the matrix B, and the number of columns in the 00088 *> matrices X and C. N >= 0. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] A 00092 *> \verbatim 00093 *> A is COMPLEX*16 array, dimension (LDA,M) 00094 *> The upper triangular matrix A. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] LDA 00098 *> \verbatim 00099 *> LDA is INTEGER 00100 *> The leading dimension of the array A. LDA >= max(1,M). 00101 *> \endverbatim 00102 *> 00103 *> \param[in] B 00104 *> \verbatim 00105 *> B is COMPLEX*16 array, dimension (LDB,N) 00106 *> The upper triangular matrix B. 00107 *> \endverbatim 00108 *> 00109 *> \param[in] LDB 00110 *> \verbatim 00111 *> LDB is INTEGER 00112 *> The leading dimension of the array B. LDB >= max(1,N). 00113 *> \endverbatim 00114 *> 00115 *> \param[in,out] C 00116 *> \verbatim 00117 *> C is COMPLEX*16 array, dimension (LDC,N) 00118 *> On entry, the M-by-N right hand side matrix C. 00119 *> On exit, C is overwritten by the solution matrix X. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] LDC 00123 *> \verbatim 00124 *> LDC is INTEGER 00125 *> The leading dimension of the array C. LDC >= max(1,M) 00126 *> \endverbatim 00127 *> 00128 *> \param[out] SCALE 00129 *> \verbatim 00130 *> SCALE is DOUBLE PRECISION 00131 *> The scale factor, scale, set <= 1 to avoid overflow in X. 00132 *> \endverbatim 00133 *> 00134 *> \param[out] INFO 00135 *> \verbatim 00136 *> INFO is INTEGER 00137 *> = 0: successful exit 00138 *> < 0: if INFO = -i, the i-th argument had an illegal value 00139 *> = 1: A and B have common or very close eigenvalues; perturbed 00140 *> values were used to solve the equation (but the matrices 00141 *> A and B are unchanged). 00142 *> \endverbatim 00143 * 00144 * Authors: 00145 * ======== 00146 * 00147 *> \author Univ. of Tennessee 00148 *> \author Univ. of California Berkeley 00149 *> \author Univ. of Colorado Denver 00150 *> \author NAG Ltd. 00151 * 00152 *> \date November 2011 00153 * 00154 *> \ingroup complex16SYcomputational 00155 * 00156 * ===================================================================== 00157 SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 00158 $ LDC, SCALE, INFO ) 00159 * 00160 * -- LAPACK computational routine (version 3.4.0) -- 00161 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00162 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00163 * November 2011 00164 * 00165 * .. Scalar Arguments .. 00166 CHARACTER TRANA, TRANB 00167 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 00168 DOUBLE PRECISION SCALE 00169 * .. 00170 * .. Array Arguments .. 00171 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) 00172 * .. 00173 * 00174 * ===================================================================== 00175 * 00176 * .. Parameters .. 00177 DOUBLE PRECISION ONE 00178 PARAMETER ( ONE = 1.0D+0 ) 00179 * .. 00180 * .. Local Scalars .. 00181 LOGICAL NOTRNA, NOTRNB 00182 INTEGER J, K, L 00183 DOUBLE PRECISION BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 00184 $ SMLNUM 00185 COMPLEX*16 A11, SUML, SUMR, VEC, X11 00186 * .. 00187 * .. Local Arrays .. 00188 DOUBLE PRECISION DUM( 1 ) 00189 * .. 00190 * .. External Functions .. 00191 LOGICAL LSAME 00192 DOUBLE PRECISION DLAMCH, ZLANGE 00193 COMPLEX*16 ZDOTC, ZDOTU, ZLADIV 00194 EXTERNAL LSAME, DLAMCH, ZLANGE, ZDOTC, ZDOTU, ZLADIV 00195 * .. 00196 * .. External Subroutines .. 00197 EXTERNAL DLABAD, XERBLA, ZDSCAL 00198 * .. 00199 * .. Intrinsic Functions .. 00200 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN 00201 * .. 00202 * .. Executable Statements .. 00203 * 00204 * Decode and Test input parameters 00205 * 00206 NOTRNA = LSAME( TRANA, 'N' ) 00207 NOTRNB = LSAME( TRANB, 'N' ) 00208 * 00209 INFO = 0 00210 IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'C' ) ) THEN 00211 INFO = -1 00212 ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'C' ) ) THEN 00213 INFO = -2 00214 ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN 00215 INFO = -3 00216 ELSE IF( M.LT.0 ) THEN 00217 INFO = -4 00218 ELSE IF( N.LT.0 ) THEN 00219 INFO = -5 00220 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00221 INFO = -7 00222 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00223 INFO = -9 00224 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00225 INFO = -11 00226 END IF 00227 IF( INFO.NE.0 ) THEN 00228 CALL XERBLA( 'ZTRSYL', -INFO ) 00229 RETURN 00230 END IF 00231 * 00232 * Quick return if possible 00233 * 00234 SCALE = ONE 00235 IF( M.EQ.0 .OR. N.EQ.0 ) 00236 $ RETURN 00237 * 00238 * Set constants to control overflow 00239 * 00240 EPS = DLAMCH( 'P' ) 00241 SMLNUM = DLAMCH( 'S' ) 00242 BIGNUM = ONE / SMLNUM 00243 CALL DLABAD( SMLNUM, BIGNUM ) 00244 SMLNUM = SMLNUM*DBLE( M*N ) / EPS 00245 BIGNUM = ONE / SMLNUM 00246 SMIN = MAX( SMLNUM, EPS*ZLANGE( 'M', M, M, A, LDA, DUM ), 00247 $ EPS*ZLANGE( 'M', N, N, B, LDB, DUM ) ) 00248 SGN = ISGN 00249 * 00250 IF( NOTRNA .AND. NOTRNB ) THEN 00251 * 00252 * Solve A*X + ISGN*X*B = scale*C. 00253 * 00254 * The (K,L)th block of X is determined starting from 00255 * bottom-left corner column by column by 00256 * 00257 * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 00258 * 00259 * Where 00260 * M L-1 00261 * R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]. 00262 * I=K+1 J=1 00263 * 00264 DO 30 L = 1, N 00265 DO 20 K = M, 1, -1 00266 * 00267 SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA, 00268 $ C( MIN( K+1, M ), L ), 1 ) 00269 SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 ) 00270 VEC = C( K, L ) - ( SUML+SGN*SUMR ) 00271 * 00272 SCALOC = ONE 00273 A11 = A( K, K ) + SGN*B( L, L ) 00274 DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 00275 IF( DA11.LE.SMIN ) THEN 00276 A11 = SMIN 00277 DA11 = SMIN 00278 INFO = 1 00279 END IF 00280 DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 00281 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00282 IF( DB.GT.BIGNUM*DA11 ) 00283 $ SCALOC = ONE / DB 00284 END IF 00285 X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 00286 * 00287 IF( SCALOC.NE.ONE ) THEN 00288 DO 10 J = 1, N 00289 CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 00290 10 CONTINUE 00291 SCALE = SCALE*SCALOC 00292 END IF 00293 C( K, L ) = X11 00294 * 00295 20 CONTINUE 00296 30 CONTINUE 00297 * 00298 ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN 00299 * 00300 * Solve A**H *X + ISGN*X*B = scale*C. 00301 * 00302 * The (K,L)th block of X is determined starting from 00303 * upper-left corner column by column by 00304 * 00305 * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 00306 * 00307 * Where 00308 * K-1 L-1 00309 * R(K,L) = SUM [A**H(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)] 00310 * I=1 J=1 00311 * 00312 DO 60 L = 1, N 00313 DO 50 K = 1, M 00314 * 00315 SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 ) 00316 SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 ) 00317 VEC = C( K, L ) - ( SUML+SGN*SUMR ) 00318 * 00319 SCALOC = ONE 00320 A11 = DCONJG( A( K, K ) ) + SGN*B( L, L ) 00321 DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 00322 IF( DA11.LE.SMIN ) THEN 00323 A11 = SMIN 00324 DA11 = SMIN 00325 INFO = 1 00326 END IF 00327 DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 00328 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00329 IF( DB.GT.BIGNUM*DA11 ) 00330 $ SCALOC = ONE / DB 00331 END IF 00332 * 00333 X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 00334 * 00335 IF( SCALOC.NE.ONE ) THEN 00336 DO 40 J = 1, N 00337 CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 00338 40 CONTINUE 00339 SCALE = SCALE*SCALOC 00340 END IF 00341 C( K, L ) = X11 00342 * 00343 50 CONTINUE 00344 60 CONTINUE 00345 * 00346 ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN 00347 * 00348 * Solve A**H*X + ISGN*X*B**H = C. 00349 * 00350 * The (K,L)th block of X is determined starting from 00351 * upper-right corner column by column by 00352 * 00353 * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L) 00354 * 00355 * Where 00356 * K-1 00357 * R(K,L) = SUM [A**H(I,K)*X(I,L)] + 00358 * I=1 00359 * N 00360 * ISGN*SUM [X(K,J)*B**H(L,J)]. 00361 * J=L+1 00362 * 00363 DO 90 L = N, 1, -1 00364 DO 80 K = 1, M 00365 * 00366 SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 ) 00367 SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC, 00368 $ B( L, MIN( L+1, N ) ), LDB ) 00369 VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) ) 00370 * 00371 SCALOC = ONE 00372 A11 = DCONJG( A( K, K )+SGN*B( L, L ) ) 00373 DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 00374 IF( DA11.LE.SMIN ) THEN 00375 A11 = SMIN 00376 DA11 = SMIN 00377 INFO = 1 00378 END IF 00379 DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 00380 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00381 IF( DB.GT.BIGNUM*DA11 ) 00382 $ SCALOC = ONE / DB 00383 END IF 00384 * 00385 X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 00386 * 00387 IF( SCALOC.NE.ONE ) THEN 00388 DO 70 J = 1, N 00389 CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 00390 70 CONTINUE 00391 SCALE = SCALE*SCALOC 00392 END IF 00393 C( K, L ) = X11 00394 * 00395 80 CONTINUE 00396 90 CONTINUE 00397 * 00398 ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN 00399 * 00400 * Solve A*X + ISGN*X*B**H = C. 00401 * 00402 * The (K,L)th block of X is determined starting from 00403 * bottom-left corner column by column by 00404 * 00405 * A(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L) 00406 * 00407 * Where 00408 * M N 00409 * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B**H(L,J)] 00410 * I=K+1 J=L+1 00411 * 00412 DO 120 L = N, 1, -1 00413 DO 110 K = M, 1, -1 00414 * 00415 SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA, 00416 $ C( MIN( K+1, M ), L ), 1 ) 00417 SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC, 00418 $ B( L, MIN( L+1, N ) ), LDB ) 00419 VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) ) 00420 * 00421 SCALOC = ONE 00422 A11 = A( K, K ) + SGN*DCONJG( B( L, L ) ) 00423 DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) ) 00424 IF( DA11.LE.SMIN ) THEN 00425 A11 = SMIN 00426 DA11 = SMIN 00427 INFO = 1 00428 END IF 00429 DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) ) 00430 IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 00431 IF( DB.GT.BIGNUM*DA11 ) 00432 $ SCALOC = ONE / DB 00433 END IF 00434 * 00435 X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 ) 00436 * 00437 IF( SCALOC.NE.ONE ) THEN 00438 DO 100 J = 1, N 00439 CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 ) 00440 100 CONTINUE 00441 SCALE = SCALE*SCALOC 00442 END IF 00443 C( K, L ) = X11 00444 * 00445 110 CONTINUE 00446 120 CONTINUE 00447 * 00448 END IF 00449 * 00450 RETURN 00451 * 00452 * End of ZTRSYL 00453 * 00454 END