Annotation of rpl/lapack/lapack/ztgsy2.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgsy2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack