![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztrsyl.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to ztrsyl.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
2: $ LDC, SCALE, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER TRANA, TRANB
11: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
12: DOUBLE PRECISION SCALE
13: * ..
14: * .. Array Arguments ..
15: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * ZTRSYL solves the complex Sylvester matrix equation:
22: *
23: * op(A)*X + X*op(B) = scale*C or
24: * op(A)*X - X*op(B) = scale*C,
25: *
26: * where op(A) = A or A**H, and A and B are both upper triangular. A is
27: * M-by-M and B is N-by-N; the right hand side C and the solution X are
28: * M-by-N; and scale is an output scale factor, set <= 1 to avoid
29: * overflow in X.
30: *
31: * Arguments
32: * =========
33: *
34: * TRANA (input) CHARACTER*1
35: * Specifies the option op(A):
36: * = 'N': op(A) = A (No transpose)
37: * = 'C': op(A) = A**H (Conjugate transpose)
38: *
39: * TRANB (input) CHARACTER*1
40: * Specifies the option op(B):
41: * = 'N': op(B) = B (No transpose)
42: * = 'C': op(B) = B**H (Conjugate transpose)
43: *
44: * ISGN (input) INTEGER
45: * Specifies the sign in the equation:
46: * = +1: solve op(A)*X + X*op(B) = scale*C
47: * = -1: solve op(A)*X - X*op(B) = scale*C
48: *
49: * M (input) INTEGER
50: * The order of the matrix A, and the number of rows in the
51: * matrices X and C. M >= 0.
52: *
53: * N (input) INTEGER
54: * The order of the matrix B, and the number of columns in the
55: * matrices X and C. N >= 0.
56: *
57: * A (input) COMPLEX*16 array, dimension (LDA,M)
58: * The upper triangular matrix A.
59: *
60: * LDA (input) INTEGER
61: * The leading dimension of the array A. LDA >= max(1,M).
62: *
63: * B (input) COMPLEX*16 array, dimension (LDB,N)
64: * The upper triangular matrix B.
65: *
66: * LDB (input) INTEGER
67: * The leading dimension of the array B. LDB >= max(1,N).
68: *
69: * C (input/output) COMPLEX*16 array, dimension (LDC,N)
70: * On entry, the M-by-N right hand side matrix C.
71: * On exit, C is overwritten by the solution matrix X.
72: *
73: * LDC (input) INTEGER
74: * The leading dimension of the array C. LDC >= max(1,M)
75: *
76: * SCALE (output) DOUBLE PRECISION
77: * The scale factor, scale, set <= 1 to avoid overflow in X.
78: *
79: * INFO (output) INTEGER
80: * = 0: successful exit
81: * < 0: if INFO = -i, the i-th argument had an illegal value
82: * = 1: A and B have common or very close eigenvalues; perturbed
83: * values were used to solve the equation (but the matrices
84: * A and B are unchanged).
85: *
86: * =====================================================================
87: *
88: * .. Parameters ..
89: DOUBLE PRECISION ONE
90: PARAMETER ( ONE = 1.0D+0 )
91: * ..
92: * .. Local Scalars ..
93: LOGICAL NOTRNA, NOTRNB
94: INTEGER J, K, L
95: DOUBLE PRECISION BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
96: $ SMLNUM
97: COMPLEX*16 A11, SUML, SUMR, VEC, X11
98: * ..
99: * .. Local Arrays ..
100: DOUBLE PRECISION DUM( 1 )
101: * ..
102: * .. External Functions ..
103: LOGICAL LSAME
104: DOUBLE PRECISION DLAMCH, ZLANGE
105: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
106: EXTERNAL LSAME, DLAMCH, ZLANGE, ZDOTC, ZDOTU, ZLADIV
107: * ..
108: * .. External Subroutines ..
109: EXTERNAL DLABAD, XERBLA, ZDSCAL
110: * ..
111: * .. Intrinsic Functions ..
112: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
113: * ..
114: * .. Executable Statements ..
115: *
116: * Decode and Test input parameters
117: *
118: NOTRNA = LSAME( TRANA, 'N' )
119: NOTRNB = LSAME( TRANB, 'N' )
120: *
121: INFO = 0
122: IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'C' ) ) THEN
123: INFO = -1
124: ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'C' ) ) THEN
125: INFO = -2
126: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
127: INFO = -3
128: ELSE IF( M.LT.0 ) THEN
129: INFO = -4
130: ELSE IF( N.LT.0 ) THEN
131: INFO = -5
132: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
133: INFO = -7
134: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
135: INFO = -9
136: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
137: INFO = -11
138: END IF
139: IF( INFO.NE.0 ) THEN
140: CALL XERBLA( 'ZTRSYL', -INFO )
141: RETURN
142: END IF
143: *
144: * Quick return if possible
145: *
146: SCALE = ONE
147: IF( M.EQ.0 .OR. N.EQ.0 )
148: $ RETURN
149: *
150: * Set constants to control overflow
151: *
152: EPS = DLAMCH( 'P' )
153: SMLNUM = DLAMCH( 'S' )
154: BIGNUM = ONE / SMLNUM
155: CALL DLABAD( SMLNUM, BIGNUM )
156: SMLNUM = SMLNUM*DBLE( M*N ) / EPS
157: BIGNUM = ONE / SMLNUM
158: SMIN = MAX( SMLNUM, EPS*ZLANGE( 'M', M, M, A, LDA, DUM ),
159: $ EPS*ZLANGE( 'M', N, N, B, LDB, DUM ) )
160: SGN = ISGN
161: *
162: IF( NOTRNA .AND. NOTRNB ) THEN
163: *
164: * Solve A*X + ISGN*X*B = scale*C.
165: *
166: * The (K,L)th block of X is determined starting from
167: * bottom-left corner column by column by
168: *
169: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
170: *
171: * Where
172: * M L-1
173: * R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)].
174: * I=K+1 J=1
175: *
176: DO 30 L = 1, N
177: DO 20 K = M, 1, -1
178: *
179: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
180: $ C( MIN( K+1, M ), L ), 1 )
181: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
182: VEC = C( K, L ) - ( SUML+SGN*SUMR )
183: *
184: SCALOC = ONE
185: A11 = A( K, K ) + SGN*B( L, L )
186: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
187: IF( DA11.LE.SMIN ) THEN
188: A11 = SMIN
189: DA11 = SMIN
190: INFO = 1
191: END IF
192: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
193: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
194: IF( DB.GT.BIGNUM*DA11 )
195: $ SCALOC = ONE / DB
196: END IF
197: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
198: *
199: IF( SCALOC.NE.ONE ) THEN
200: DO 10 J = 1, N
201: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
202: 10 CONTINUE
203: SCALE = SCALE*SCALOC
204: END IF
205: C( K, L ) = X11
206: *
207: 20 CONTINUE
208: 30 CONTINUE
209: *
210: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
211: *
212: * Solve A' *X + ISGN*X*B = scale*C.
213: *
214: * The (K,L)th block of X is determined starting from
215: * upper-left corner column by column by
216: *
217: * A'(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
218: *
219: * Where
220: * K-1 L-1
221: * R(K,L) = SUM [A'(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]
222: * I=1 J=1
223: *
224: DO 60 L = 1, N
225: DO 50 K = 1, M
226: *
227: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
228: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
229: VEC = C( K, L ) - ( SUML+SGN*SUMR )
230: *
231: SCALOC = ONE
232: A11 = DCONJG( A( K, K ) ) + SGN*B( L, L )
233: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
234: IF( DA11.LE.SMIN ) THEN
235: A11 = SMIN
236: DA11 = SMIN
237: INFO = 1
238: END IF
239: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
240: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
241: IF( DB.GT.BIGNUM*DA11 )
242: $ SCALOC = ONE / DB
243: END IF
244: *
245: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
246: *
247: IF( SCALOC.NE.ONE ) THEN
248: DO 40 J = 1, N
249: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
250: 40 CONTINUE
251: SCALE = SCALE*SCALOC
252: END IF
253: C( K, L ) = X11
254: *
255: 50 CONTINUE
256: 60 CONTINUE
257: *
258: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
259: *
260: * Solve A'*X + ISGN*X*B' = C.
261: *
262: * The (K,L)th block of X is determined starting from
263: * upper-right corner column by column by
264: *
265: * A'(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
266: *
267: * Where
268: * K-1
269: * R(K,L) = SUM [A'(I,K)*X(I,L)] +
270: * I=1
271: * N
272: * ISGN*SUM [X(K,J)*B'(L,J)].
273: * J=L+1
274: *
275: DO 90 L = N, 1, -1
276: DO 80 K = 1, M
277: *
278: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
279: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
280: $ B( L, MIN( L+1, N ) ), LDB )
281: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
282: *
283: SCALOC = ONE
284: A11 = DCONJG( A( K, K )+SGN*B( L, L ) )
285: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
286: IF( DA11.LE.SMIN ) THEN
287: A11 = SMIN
288: DA11 = SMIN
289: INFO = 1
290: END IF
291: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
292: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
293: IF( DB.GT.BIGNUM*DA11 )
294: $ SCALOC = ONE / DB
295: END IF
296: *
297: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
298: *
299: IF( SCALOC.NE.ONE ) THEN
300: DO 70 J = 1, N
301: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
302: 70 CONTINUE
303: SCALE = SCALE*SCALOC
304: END IF
305: C( K, L ) = X11
306: *
307: 80 CONTINUE
308: 90 CONTINUE
309: *
310: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
311: *
312: * Solve A*X + ISGN*X*B' = C.
313: *
314: * The (K,L)th block of X is determined starting from
315: * bottom-left corner column by column by
316: *
317: * A(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
318: *
319: * Where
320: * M N
321: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B'(L,J)]
322: * I=K+1 J=L+1
323: *
324: DO 120 L = N, 1, -1
325: DO 110 K = M, 1, -1
326: *
327: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
328: $ C( MIN( K+1, M ), L ), 1 )
329: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
330: $ B( L, MIN( L+1, N ) ), LDB )
331: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
332: *
333: SCALOC = ONE
334: A11 = A( K, K ) + SGN*DCONJG( B( L, L ) )
335: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
336: IF( DA11.LE.SMIN ) THEN
337: A11 = SMIN
338: DA11 = SMIN
339: INFO = 1
340: END IF
341: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
342: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
343: IF( DB.GT.BIGNUM*DA11 )
344: $ SCALOC = ONE / DB
345: END IF
346: *
347: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
348: *
349: IF( SCALOC.NE.ONE ) THEN
350: DO 100 J = 1, N
351: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
352: 100 CONTINUE
353: SCALE = SCALE*SCALOC
354: END IF
355: C( K, L ) = X11
356: *
357: 110 CONTINUE
358: 120 CONTINUE
359: *
360: END IF
361: *
362: RETURN
363: *
364: * End of ZTRSYL
365: *
366: END