![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTGSY2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTGSY2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsy2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsy2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsy2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 00022 * LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 00023 * INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER TRANS 00027 * INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N 00028 * DOUBLE PRECISION RDSCAL, RDSUM, SCALE 00029 * .. 00030 * .. Array Arguments .. 00031 * COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), 00032 * $ D( LDD, * ), E( LDE, * ), F( LDF, * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> ZTGSY2 solves the generalized Sylvester equation 00042 *> 00043 *> A * R - L * B = scale * C (1) 00044 *> D * R - L * E = scale * F 00045 *> 00046 *> using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices, 00047 *> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, 00048 *> N-by-N and M-by-N, respectively. A, B, D and E are upper triangular 00049 *> (i.e., (A,D) and (B,E) in generalized Schur form). 00050 *> 00051 *> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output 00052 *> scaling factor chosen to avoid overflow. 00053 *> 00054 *> In matrix notation solving equation (1) corresponds to solve 00055 *> Zx = scale * b, where Z is defined as 00056 *> 00057 *> Z = [ kron(In, A) -kron(B**H, Im) ] (2) 00058 *> [ kron(In, D) -kron(E**H, Im) ], 00059 *> 00060 *> Ik is the identity matrix of size k and X**H is the conjuguate transpose of X. 00061 *> kron(X, Y) is the Kronecker product between the matrices X and Y. 00062 *> 00063 *> If TRANS = 'C', y in the conjugate transposed system Z**H*y = scale*b 00064 *> is solved for, which is equivalent to solve for R and L in 00065 *> 00066 *> A**H * R + D**H * L = scale * C (3) 00067 *> R * B**H + L * E**H = scale * -F 00068 *> 00069 *> This case is used to compute an estimate of Dif[(A, D), (B, E)] = 00070 *> = sigma_min(Z) using reverse communicaton with ZLACON. 00071 *> 00072 *> ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL 00073 *> of an upper bound on the separation between to matrix pairs. Then 00074 *> the input (A, D), (B, E) are sub-pencils of two matrix pairs in 00075 *> ZTGSYL. 00076 *> \endverbatim 00077 * 00078 * Arguments: 00079 * ========== 00080 * 00081 *> \param[in] TRANS 00082 *> \verbatim 00083 *> TRANS is CHARACTER*1 00084 *> = 'N', solve the generalized Sylvester equation (1). 00085 *> = 'T': solve the 'transposed' system (3). 00086 *> \endverbatim 00087 *> 00088 *> \param[in] IJOB 00089 *> \verbatim 00090 *> IJOB is INTEGER 00091 *> Specifies what kind of functionality to be performed. 00092 *> =0: solve (1) only. 00093 *> =1: A contribution from this subsystem to a Frobenius 00094 *> norm-based estimate of the separation between two matrix 00095 *> pairs is computed. (look ahead strategy is used). 00096 *> =2: A contribution from this subsystem to a Frobenius 00097 *> norm-based estimate of the separation between two matrix 00098 *> pairs is computed. (DGECON on sub-systems is used.) 00099 *> Not referenced if TRANS = 'T'. 00100 *> \endverbatim 00101 *> 00102 *> \param[in] M 00103 *> \verbatim 00104 *> M is INTEGER 00105 *> On entry, M specifies the order of A and D, and the row 00106 *> dimension of C, F, R and L. 00107 *> \endverbatim 00108 *> 00109 *> \param[in] N 00110 *> \verbatim 00111 *> N is INTEGER 00112 *> On entry, N specifies the order of B and E, and the column 00113 *> dimension of C, F, R and L. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] A 00117 *> \verbatim 00118 *> A is COMPLEX*16 array, dimension (LDA, M) 00119 *> On entry, A contains an upper triangular matrix. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] LDA 00123 *> \verbatim 00124 *> LDA is INTEGER 00125 *> The leading dimension of the matrix A. LDA >= max(1, M). 00126 *> \endverbatim 00127 *> 00128 *> \param[in] B 00129 *> \verbatim 00130 *> B is COMPLEX*16 array, dimension (LDB, N) 00131 *> On entry, B contains an upper triangular matrix. 00132 *> \endverbatim 00133 *> 00134 *> \param[in] LDB 00135 *> \verbatim 00136 *> LDB is INTEGER 00137 *> The leading dimension of the matrix B. LDB >= max(1, N). 00138 *> \endverbatim 00139 *> 00140 *> \param[in,out] C 00141 *> \verbatim 00142 *> C is COMPLEX*16 array, dimension (LDC, N) 00143 *> On entry, C contains the right-hand-side of the first matrix 00144 *> equation in (1). 00145 *> On exit, if IJOB = 0, C has been overwritten by the solution 00146 *> R. 00147 *> \endverbatim 00148 *> 00149 *> \param[in] LDC 00150 *> \verbatim 00151 *> LDC is INTEGER 00152 *> The leading dimension of the matrix C. LDC >= max(1, M). 00153 *> \endverbatim 00154 *> 00155 *> \param[in] D 00156 *> \verbatim 00157 *> D is COMPLEX*16 array, dimension (LDD, M) 00158 *> On entry, D contains an upper triangular matrix. 00159 *> \endverbatim 00160 *> 00161 *> \param[in] LDD 00162 *> \verbatim 00163 *> LDD is INTEGER 00164 *> The leading dimension of the matrix D. LDD >= max(1, M). 00165 *> \endverbatim 00166 *> 00167 *> \param[in] E 00168 *> \verbatim 00169 *> E is COMPLEX*16 array, dimension (LDE, N) 00170 *> On entry, E contains an upper triangular matrix. 00171 *> \endverbatim 00172 *> 00173 *> \param[in] LDE 00174 *> \verbatim 00175 *> LDE is INTEGER 00176 *> The leading dimension of the matrix E. LDE >= max(1, N). 00177 *> \endverbatim 00178 *> 00179 *> \param[in,out] F 00180 *> \verbatim 00181 *> F is COMPLEX*16 array, dimension (LDF, N) 00182 *> On entry, F contains the right-hand-side of the second matrix 00183 *> equation in (1). 00184 *> On exit, if IJOB = 0, F has been overwritten by the solution 00185 *> L. 00186 *> \endverbatim 00187 *> 00188 *> \param[in] LDF 00189 *> \verbatim 00190 *> LDF is INTEGER 00191 *> The leading dimension of the matrix F. LDF >= max(1, M). 00192 *> \endverbatim 00193 *> 00194 *> \param[out] SCALE 00195 *> \verbatim 00196 *> SCALE is DOUBLE PRECISION 00197 *> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions 00198 *> R and L (C and F on entry) will hold the solutions to a 00199 *> slightly perturbed system but the input matrices A, B, D and 00200 *> E have not been changed. If SCALE = 0, R and L will hold the 00201 *> solutions to the homogeneous system with C = F = 0. 00202 *> Normally, SCALE = 1. 00203 *> \endverbatim 00204 *> 00205 *> \param[in,out] RDSUM 00206 *> \verbatim 00207 *> RDSUM is DOUBLE PRECISION 00208 *> On entry, the sum of squares of computed contributions to 00209 *> the Dif-estimate under computation by ZTGSYL, where the 00210 *> scaling factor RDSCAL (see below) has been factored out. 00211 *> On exit, the corresponding sum of squares updated with the 00212 *> contributions from the current sub-system. 00213 *> If TRANS = 'T' RDSUM is not touched. 00214 *> NOTE: RDSUM only makes sense when ZTGSY2 is called by 00215 *> ZTGSYL. 00216 *> \endverbatim 00217 *> 00218 *> \param[in,out] RDSCAL 00219 *> \verbatim 00220 *> RDSCAL is DOUBLE PRECISION 00221 *> On entry, scaling factor used to prevent overflow in RDSUM. 00222 *> On exit, RDSCAL is updated w.r.t. the current contributions 00223 *> in RDSUM. 00224 *> If TRANS = 'T', RDSCAL is not touched. 00225 *> NOTE: RDSCAL only makes sense when ZTGSY2 is called by 00226 *> ZTGSYL. 00227 *> \endverbatim 00228 *> 00229 *> \param[out] INFO 00230 *> \verbatim 00231 *> INFO is INTEGER 00232 *> On exit, if INFO is set to 00233 *> =0: Successful exit 00234 *> <0: If INFO = -i, input argument number i is illegal. 00235 *> >0: The matrix pairs (A, D) and (B, E) have common or very 00236 *> close eigenvalues. 00237 *> \endverbatim 00238 * 00239 * Authors: 00240 * ======== 00241 * 00242 *> \author Univ. of Tennessee 00243 *> \author Univ. of California Berkeley 00244 *> \author Univ. of Colorado Denver 00245 *> \author NAG Ltd. 00246 * 00247 *> \date November 2011 00248 * 00249 *> \ingroup complex16SYauxiliary 00250 * 00251 *> \par Contributors: 00252 * ================== 00253 *> 00254 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00255 *> Umea University, S-901 87 Umea, Sweden. 00256 * 00257 * ===================================================================== 00258 SUBROUTINE ZTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 00259 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 00260 $ INFO ) 00261 * 00262 * -- LAPACK auxiliary routine (version 3.4.0) -- 00263 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00264 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00265 * November 2011 00266 * 00267 * .. Scalar Arguments .. 00268 CHARACTER TRANS 00269 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N 00270 DOUBLE PRECISION RDSCAL, RDSUM, SCALE 00271 * .. 00272 * .. Array Arguments .. 00273 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), 00274 $ D( LDD, * ), E( LDE, * ), F( LDF, * ) 00275 * .. 00276 * 00277 * ===================================================================== 00278 * 00279 * .. Parameters .. 00280 DOUBLE PRECISION ZERO, ONE 00281 INTEGER LDZ 00282 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, LDZ = 2 ) 00283 * .. 00284 * .. Local Scalars .. 00285 LOGICAL NOTRAN 00286 INTEGER I, IERR, J, K 00287 DOUBLE PRECISION SCALOC 00288 COMPLEX*16 ALPHA 00289 * .. 00290 * .. Local Arrays .. 00291 INTEGER IPIV( LDZ ), JPIV( LDZ ) 00292 COMPLEX*16 RHS( LDZ ), Z( LDZ, LDZ ) 00293 * .. 00294 * .. External Functions .. 00295 LOGICAL LSAME 00296 EXTERNAL LSAME 00297 * .. 00298 * .. External Subroutines .. 00299 EXTERNAL XERBLA, ZAXPY, ZGESC2, ZGETC2, ZLATDF, ZSCAL 00300 * .. 00301 * .. Intrinsic Functions .. 00302 INTRINSIC DCMPLX, DCONJG, MAX 00303 * .. 00304 * .. Executable Statements .. 00305 * 00306 * Decode and test input parameters 00307 * 00308 INFO = 0 00309 IERR = 0 00310 NOTRAN = LSAME( TRANS, 'N' ) 00311 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00312 INFO = -1 00313 ELSE IF( NOTRAN ) THEN 00314 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN 00315 INFO = -2 00316 END IF 00317 END IF 00318 IF( INFO.EQ.0 ) THEN 00319 IF( M.LE.0 ) THEN 00320 INFO = -3 00321 ELSE IF( N.LE.0 ) THEN 00322 INFO = -4 00323 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00324 INFO = -5 00325 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00326 INFO = -8 00327 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00328 INFO = -10 00329 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 00330 INFO = -12 00331 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 00332 INFO = -14 00333 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 00334 INFO = -16 00335 END IF 00336 END IF 00337 IF( INFO.NE.0 ) THEN 00338 CALL XERBLA( 'ZTGSY2', -INFO ) 00339 RETURN 00340 END IF 00341 * 00342 IF( NOTRAN ) THEN 00343 * 00344 * Solve (I, J) - system 00345 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 00346 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 00347 * for I = M, M - 1, ..., 1; J = 1, 2, ..., N 00348 * 00349 SCALE = ONE 00350 SCALOC = ONE 00351 DO 30 J = 1, N 00352 DO 20 I = M, 1, -1 00353 * 00354 * Build 2 by 2 system 00355 * 00356 Z( 1, 1 ) = A( I, I ) 00357 Z( 2, 1 ) = D( I, I ) 00358 Z( 1, 2 ) = -B( J, J ) 00359 Z( 2, 2 ) = -E( J, J ) 00360 * 00361 * Set up right hand side(s) 00362 * 00363 RHS( 1 ) = C( I, J ) 00364 RHS( 2 ) = F( I, J ) 00365 * 00366 * Solve Z * x = RHS 00367 * 00368 CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR ) 00369 IF( IERR.GT.0 ) 00370 $ INFO = IERR 00371 IF( IJOB.EQ.0 ) THEN 00372 CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00373 IF( SCALOC.NE.ONE ) THEN 00374 DO 10 K = 1, N 00375 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 00376 $ C( 1, K ), 1 ) 00377 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 00378 $ F( 1, K ), 1 ) 00379 10 CONTINUE 00380 SCALE = SCALE*SCALOC 00381 END IF 00382 ELSE 00383 CALL ZLATDF( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL, 00384 $ IPIV, JPIV ) 00385 END IF 00386 * 00387 * Unpack solution vector(s) 00388 * 00389 C( I, J ) = RHS( 1 ) 00390 F( I, J ) = RHS( 2 ) 00391 * 00392 * Substitute R(I, J) and L(I, J) into remaining equation. 00393 * 00394 IF( I.GT.1 ) THEN 00395 ALPHA = -RHS( 1 ) 00396 CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 ) 00397 CALL ZAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 ) 00398 END IF 00399 IF( J.LT.N ) THEN 00400 CALL ZAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB, 00401 $ C( I, J+1 ), LDC ) 00402 CALL ZAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE, 00403 $ F( I, J+1 ), LDF ) 00404 END IF 00405 * 00406 20 CONTINUE 00407 30 CONTINUE 00408 ELSE 00409 * 00410 * Solve transposed (I, J) - system: 00411 * A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J) 00412 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 00413 * for I = 1, 2, ..., M, J = N, N - 1, ..., 1 00414 * 00415 SCALE = ONE 00416 SCALOC = ONE 00417 DO 80 I = 1, M 00418 DO 70 J = N, 1, -1 00419 * 00420 * Build 2 by 2 system Z**H 00421 * 00422 Z( 1, 1 ) = DCONJG( A( I, I ) ) 00423 Z( 2, 1 ) = -DCONJG( B( J, J ) ) 00424 Z( 1, 2 ) = DCONJG( D( I, I ) ) 00425 Z( 2, 2 ) = -DCONJG( E( J, J ) ) 00426 * 00427 * 00428 * Set up right hand side(s) 00429 * 00430 RHS( 1 ) = C( I, J ) 00431 RHS( 2 ) = F( I, J ) 00432 * 00433 * Solve Z**H * x = RHS 00434 * 00435 CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR ) 00436 IF( IERR.GT.0 ) 00437 $ INFO = IERR 00438 CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00439 IF( SCALOC.NE.ONE ) THEN 00440 DO 40 K = 1, N 00441 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ), 00442 $ 1 ) 00443 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ), 00444 $ 1 ) 00445 40 CONTINUE 00446 SCALE = SCALE*SCALOC 00447 END IF 00448 * 00449 * Unpack solution vector(s) 00450 * 00451 C( I, J ) = RHS( 1 ) 00452 F( I, J ) = RHS( 2 ) 00453 * 00454 * Substitute R(I, J) and L(I, J) into remaining equation. 00455 * 00456 DO 50 K = 1, J - 1 00457 F( I, K ) = F( I, K ) + RHS( 1 )*DCONJG( B( K, J ) ) + 00458 $ RHS( 2 )*DCONJG( E( K, J ) ) 00459 50 CONTINUE 00460 DO 60 K = I + 1, M 00461 C( K, J ) = C( K, J ) - DCONJG( A( I, K ) )*RHS( 1 ) - 00462 $ DCONJG( D( I, K ) )*RHS( 2 ) 00463 60 CONTINUE 00464 * 00465 70 CONTINUE 00466 80 CONTINUE 00467 END IF 00468 RETURN 00469 * 00470 * End of ZTGSY2 00471 * 00472 END