1 SUBROUTINE ZTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 2 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 3 $ INFO ) 4* 5* -- LAPACK auxiliary routine (version 3.2) -- 6* -- LAPACK is a software package provided by Univ. of Tennessee, -- 7* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 8* November 2006 9* 10* .. Scalar Arguments .. 11 CHARACTER TRANS 12 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N 13 DOUBLE PRECISION RDSCAL, RDSUM, SCALE 14* .. 15* .. Array Arguments .. 16 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), 17 $ D( LDD, * ), E( LDE, * ), F( LDF, * ) 18* .. 19* 20* Purpose 21* ======= 22* 23* ZTGSY2 solves the generalized Sylvester equation 24* 25* A * R - L * B = scale * C (1) 26* D * R - L * E = scale * F 27* 28* using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices, 29* (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, 30* N-by-N and M-by-N, respectively. A, B, D and E are upper triangular 31* (i.e., (A,D) and (B,E) in generalized Schur form). 32* 33* The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output 34* scaling factor chosen to avoid overflow. 35* 36* In matrix notation solving equation (1) corresponds to solve 37* Zx = scale * b, where Z is defined as 38* 39* Z = [ kron(In, A) -kron(B', Im) ] (2) 40* [ kron(In, D) -kron(E', Im) ], 41* 42* Ik is the identity matrix of size k and X' is the transpose of X. 43* kron(X, Y) is the Kronecker product between the matrices X and Y. 44* 45* If TRANS = 'C', y in the conjugate transposed system Z'y = scale*b 46* is solved for, which is equivalent to solve for R and L in 47* 48* A' * R + D' * L = scale * C (3) 49* R * B' + L * E' = scale * -F 50* 51* This case is used to compute an estimate of Dif[(A, D), (B, E)] = 52* = sigma_min(Z) using reverse communicaton with ZLACON. 53* 54* ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL 55* of an upper bound on the separation between to matrix pairs. Then 56* the input (A, D), (B, E) are sub-pencils of two matrix pairs in 57* ZTGSYL. 58* 59* Arguments 60* ========= 61* 62* TRANS (input) CHARACTER*1 63* = 'N', solve the generalized Sylvester equation (1). 64* = 'T': solve the 'transposed' system (3). 65* 66* IJOB (input) INTEGER 67* Specifies what kind of functionality to be performed. 68* =0: solve (1) only. 69* =1: A contribution from this subsystem to a Frobenius 70* norm-based estimate of the separation between two matrix 71* pairs is computed. (look ahead strategy is used). 72* =2: A contribution from this subsystem to a Frobenius 73* norm-based estimate of the separation between two matrix 74* pairs is computed. (DGECON on sub-systems is used.) 75* Not referenced if TRANS = 'T'. 76* 77* M (input) INTEGER 78* On entry, M specifies the order of A and D, and the row 79* dimension of C, F, R and L. 80* 81* N (input) INTEGER 82* On entry, N specifies the order of B and E, and the column 83* dimension of C, F, R and L. 84* 85* A (input) COMPLEX*16 array, dimension (LDA, M) 86* On entry, A contains an upper triangular matrix. 87* 88* LDA (input) INTEGER 89* The leading dimension of the matrix A. LDA >= max(1, M). 90* 91* B (input) COMPLEX*16 array, dimension (LDB, N) 92* On entry, B contains an upper triangular matrix. 93* 94* LDB (input) INTEGER 95* The leading dimension of the matrix B. LDB >= max(1, N). 96* 97* C (input/output) COMPLEX*16 array, dimension (LDC, N) 98* On entry, C contains the right-hand-side of the first matrix 99* equation in (1). 100* On exit, if IJOB = 0, C has been overwritten by the solution 101* R. 102* 103* LDC (input) INTEGER 104* The leading dimension of the matrix C. LDC >= max(1, M). 105* 106* D (input) COMPLEX*16 array, dimension (LDD, M) 107* On entry, D contains an upper triangular matrix. 108* 109* LDD (input) INTEGER 110* The leading dimension of the matrix D. LDD >= max(1, M). 111* 112* E (input) COMPLEX*16 array, dimension (LDE, N) 113* On entry, E contains an upper triangular matrix. 114* 115* LDE (input) INTEGER 116* The leading dimension of the matrix E. LDE >= max(1, N). 117* 118* F (input/output) COMPLEX*16 array, dimension (LDF, N) 119* On entry, F contains the right-hand-side of the second matrix 120* equation in (1). 121* On exit, if IJOB = 0, F has been overwritten by the solution 122* L. 123* 124* LDF (input) INTEGER 125* The leading dimension of the matrix F. LDF >= max(1, M). 126* 127* SCALE (output) DOUBLE PRECISION 128* On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions 129* R and L (C and F on entry) will hold the solutions to a 130* slightly perturbed system but the input matrices A, B, D and 131* E have not been changed. If SCALE = 0, R and L will hold the 132* solutions to the homogeneous system with C = F = 0. 133* Normally, SCALE = 1. 134* 135* RDSUM (input/output) DOUBLE PRECISION 136* On entry, the sum of squares of computed contributions to 137* the Dif-estimate under computation by ZTGSYL, where the 138* scaling factor RDSCAL (see below) has been factored out. 139* On exit, the corresponding sum of squares updated with the 140* contributions from the current sub-system. 141* If TRANS = 'T' RDSUM is not touched. 142* NOTE: RDSUM only makes sense when ZTGSY2 is called by 143* ZTGSYL. 144* 145* RDSCAL (input/output) DOUBLE PRECISION 146* On entry, scaling factor used to prevent overflow in RDSUM. 147* On exit, RDSCAL is updated w.r.t. the current contributions 148* in RDSUM. 149* If TRANS = 'T', RDSCAL is not touched. 150* NOTE: RDSCAL only makes sense when ZTGSY2 is called by 151* ZTGSYL. 152* 153* INFO (output) INTEGER 154* On exit, if INFO is set to 155* =0: Successful exit 156* <0: If INFO = -i, input argument number i is illegal. 157* >0: The matrix pairs (A, D) and (B, E) have common or very 158* close eigenvalues. 159* 160* Further Details 161* =============== 162* 163* Based on contributions by 164* Bo Kagstrom and Peter Poromaa, Department of Computing Science, 165* Umea University, S-901 87 Umea, Sweden. 166* 167* ===================================================================== 168* 169* .. Parameters .. 170 DOUBLE PRECISION ZERO, ONE 171 INTEGER LDZ 172 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, LDZ = 2 ) 173* .. 174* .. Local Scalars .. 175 LOGICAL NOTRAN 176 INTEGER I, IERR, J, K 177 DOUBLE PRECISION SCALOC 178 COMPLEX*16 ALPHA 179* .. 180* .. Local Arrays .. 181 INTEGER IPIV( LDZ ), JPIV( LDZ ) 182 COMPLEX*16 RHS( LDZ ), Z( LDZ, LDZ ) 183* .. 184* .. External Functions .. 185 LOGICAL LSAME 186 EXTERNAL LSAME 187* .. 188* .. External Subroutines .. 189 EXTERNAL XERBLA, ZAXPY, ZGESC2, ZGETC2, ZLATDF, ZSCAL 190* .. 191* .. Intrinsic Functions .. 192 INTRINSIC DCMPLX, DCONJG, MAX 193* .. 194* .. Executable Statements .. 195* 196* Decode and test input parameters 197* 198 INFO = 0 199 IERR = 0 200 NOTRAN = LSAME( TRANS, 'N' ) 201 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 202 INFO = -1 203 ELSE IF( NOTRAN ) THEN 204 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN 205 INFO = -2 206 END IF 207 END IF 208 IF( INFO.EQ.0 ) THEN 209 IF( M.LE.0 ) THEN 210 INFO = -3 211 ELSE IF( N.LE.0 ) THEN 212 INFO = -4 213 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 214 INFO = -5 215 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 216 INFO = -8 217 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 218 INFO = -10 219 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 220 INFO = -12 221 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 222 INFO = -14 223 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 224 INFO = -16 225 END IF 226 END IF 227 IF( INFO.NE.0 ) THEN 228 CALL XERBLA( 'ZTGSY2', -INFO ) 229 RETURN 230 END IF 231* 232 IF( NOTRAN ) THEN 233* 234* Solve (I, J) - system 235* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 236* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 237* for I = M, M - 1, ..., 1; J = 1, 2, ..., N 238* 239 SCALE = ONE 240 SCALOC = ONE 241 DO 30 J = 1, N 242 DO 20 I = M, 1, -1 243* 244* Build 2 by 2 system 245* 246 Z( 1, 1 ) = A( I, I ) 247 Z( 2, 1 ) = D( I, I ) 248 Z( 1, 2 ) = -B( J, J ) 249 Z( 2, 2 ) = -E( J, J ) 250* 251* Set up right hand side(s) 252* 253 RHS( 1 ) = C( I, J ) 254 RHS( 2 ) = F( I, J ) 255* 256* Solve Z * x = RHS 257* 258 CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR ) 259 IF( IERR.GT.0 ) 260 $ INFO = IERR 261 IF( IJOB.EQ.0 ) THEN 262 CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 263 IF( SCALOC.NE.ONE ) THEN 264 DO 10 K = 1, N 265 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 266 $ C( 1, K ), 1 ) 267 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 268 $ F( 1, K ), 1 ) 269 10 CONTINUE 270 SCALE = SCALE*SCALOC 271 END IF 272 ELSE 273 CALL ZLATDF( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL, 274 $ IPIV, JPIV ) 275 END IF 276* 277* Unpack solution vector(s) 278* 279 C( I, J ) = RHS( 1 ) 280 F( I, J ) = RHS( 2 ) 281* 282* Substitute R(I, J) and L(I, J) into remaining equation. 283* 284 IF( I.GT.1 ) THEN 285 ALPHA = -RHS( 1 ) 286 CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 ) 287 CALL ZAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 ) 288 END IF 289 IF( J.LT.N ) THEN 290 CALL ZAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB, 291 $ C( I, J+1 ), LDC ) 292 CALL ZAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE, 293 $ F( I, J+1 ), LDF ) 294 END IF 295* 296 20 CONTINUE 297 30 CONTINUE 298 ELSE 299* 300* Solve transposed (I, J) - system: 301* A(I, I)' * R(I, J) + D(I, I)' * L(J, J) = C(I, J) 302* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 303* for I = 1, 2, ..., M, J = N, N - 1, ..., 1 304* 305 SCALE = ONE 306 SCALOC = ONE 307 DO 80 I = 1, M 308 DO 70 J = N, 1, -1 309* 310* Build 2 by 2 system Z' 311* 312 Z( 1, 1 ) = DCONJG( A( I, I ) ) 313 Z( 2, 1 ) = -DCONJG( B( J, J ) ) 314 Z( 1, 2 ) = DCONJG( D( I, I ) ) 315 Z( 2, 2 ) = -DCONJG( E( J, J ) ) 316* 317* 318* Set up right hand side(s) 319* 320 RHS( 1 ) = C( I, J ) 321 RHS( 2 ) = F( I, J ) 322* 323* Solve Z' * x = RHS 324* 325 CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR ) 326 IF( IERR.GT.0 ) 327 $ INFO = IERR 328 CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 329 IF( SCALOC.NE.ONE ) THEN 330 DO 40 K = 1, N 331 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ), 332 $ 1 ) 333 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ), 334 $ 1 ) 335 40 CONTINUE 336 SCALE = SCALE*SCALOC 337 END IF 338* 339* Unpack solution vector(s) 340* 341 C( I, J ) = RHS( 1 ) 342 F( I, J ) = RHS( 2 ) 343* 344* Substitute R(I, J) and L(I, J) into remaining equation. 345* 346 DO 50 K = 1, J - 1 347 F( I, K ) = F( I, K ) + RHS( 1 )*DCONJG( B( K, J ) ) + 348 $ RHS( 2 )*DCONJG( E( K, J ) ) 349 50 CONTINUE 350 DO 60 K = I + 1, M 351 C( K, J ) = C( K, J ) - DCONJG( A( I, K ) )*RHS( 1 ) - 352 $ DCONJG( D( I, K ) )*RHS( 2 ) 353 60 CONTINUE 354* 355 70 CONTINUE 356 80 CONTINUE 357 END IF 358 RETURN 359* 360* End of ZTGSY2 361* 362 END 363