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