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