1: *> \brief \b ZTRSYL
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTRSYL + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsyl.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsyl.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsyl.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
22: * LDC, SCALE, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER TRANA, TRANB
26: * INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
27: * DOUBLE PRECISION SCALE
28: * ..
29: * .. Array Arguments ..
30: * COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZTRSYL solves the complex Sylvester matrix equation:
40: *>
41: *> op(A)*X + X*op(B) = scale*C or
42: *> op(A)*X - X*op(B) = scale*C,
43: *>
44: *> where op(A) = A or A**H, and A and B are both upper triangular. A is
45: *> M-by-M and B is N-by-N; the right hand side C and the solution X are
46: *> M-by-N; and scale is an output scale factor, set <= 1 to avoid
47: *> overflow in X.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] TRANA
54: *> \verbatim
55: *> TRANA is CHARACTER*1
56: *> Specifies the option op(A):
57: *> = 'N': op(A) = A (No transpose)
58: *> = 'C': op(A) = A**H (Conjugate transpose)
59: *> \endverbatim
60: *>
61: *> \param[in] TRANB
62: *> \verbatim
63: *> TRANB is CHARACTER*1
64: *> Specifies the option op(B):
65: *> = 'N': op(B) = B (No transpose)
66: *> = 'C': op(B) = B**H (Conjugate transpose)
67: *> \endverbatim
68: *>
69: *> \param[in] ISGN
70: *> \verbatim
71: *> ISGN is INTEGER
72: *> Specifies the sign in the equation:
73: *> = +1: solve op(A)*X + X*op(B) = scale*C
74: *> = -1: solve op(A)*X - X*op(B) = scale*C
75: *> \endverbatim
76: *>
77: *> \param[in] M
78: *> \verbatim
79: *> M is INTEGER
80: *> The order of the matrix A, and the number of rows in the
81: *> matrices X and C. M >= 0.
82: *> \endverbatim
83: *>
84: *> \param[in] N
85: *> \verbatim
86: *> N is INTEGER
87: *> The order of the matrix B, and the number of columns in the
88: *> matrices X and C. N >= 0.
89: *> \endverbatim
90: *>
91: *> \param[in] A
92: *> \verbatim
93: *> A is COMPLEX*16 array, dimension (LDA,M)
94: *> The upper triangular matrix A.
95: *> \endverbatim
96: *>
97: *> \param[in] LDA
98: *> \verbatim
99: *> LDA is INTEGER
100: *> The leading dimension of the array A. LDA >= max(1,M).
101: *> \endverbatim
102: *>
103: *> \param[in] B
104: *> \verbatim
105: *> B is COMPLEX*16 array, dimension (LDB,N)
106: *> The upper triangular matrix B.
107: *> \endverbatim
108: *>
109: *> \param[in] LDB
110: *> \verbatim
111: *> LDB is INTEGER
112: *> The leading dimension of the array B. LDB >= max(1,N).
113: *> \endverbatim
114: *>
115: *> \param[in,out] C
116: *> \verbatim
117: *> C is COMPLEX*16 array, dimension (LDC,N)
118: *> On entry, the M-by-N right hand side matrix C.
119: *> On exit, C is overwritten by the solution matrix X.
120: *> \endverbatim
121: *>
122: *> \param[in] LDC
123: *> \verbatim
124: *> LDC is INTEGER
125: *> The leading dimension of the array C. LDC >= max(1,M)
126: *> \endverbatim
127: *>
128: *> \param[out] SCALE
129: *> \verbatim
130: *> SCALE is DOUBLE PRECISION
131: *> The scale factor, scale, set <= 1 to avoid overflow in X.
132: *> \endverbatim
133: *>
134: *> \param[out] INFO
135: *> \verbatim
136: *> INFO is INTEGER
137: *> = 0: successful exit
138: *> < 0: if INFO = -i, the i-th argument had an illegal value
139: *> = 1: A and B have common or very close eigenvalues; perturbed
140: *> values were used to solve the equation (but the matrices
141: *> A and B are unchanged).
142: *> \endverbatim
143: *
144: * Authors:
145: * ========
146: *
147: *> \author Univ. of Tennessee
148: *> \author Univ. of California Berkeley
149: *> \author Univ. of Colorado Denver
150: *> \author NAG Ltd.
151: *
152: *> \ingroup complex16SYcomputational
153: *
154: * =====================================================================
155: SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
156: $ LDC, SCALE, INFO )
157: *
158: * -- LAPACK computational routine --
159: * -- LAPACK is a software package provided by Univ. of Tennessee, --
160: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
161: *
162: * .. Scalar Arguments ..
163: CHARACTER TRANA, TRANB
164: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
165: DOUBLE PRECISION SCALE
166: * ..
167: * .. Array Arguments ..
168: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
169: * ..
170: *
171: * =====================================================================
172: *
173: * .. Parameters ..
174: DOUBLE PRECISION ONE
175: PARAMETER ( ONE = 1.0D+0 )
176: * ..
177: * .. Local Scalars ..
178: LOGICAL NOTRNA, NOTRNB
179: INTEGER J, K, L
180: DOUBLE PRECISION BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
181: $ SMLNUM
182: COMPLEX*16 A11, SUML, SUMR, VEC, X11
183: * ..
184: * .. Local Arrays ..
185: DOUBLE PRECISION DUM( 1 )
186: * ..
187: * .. External Functions ..
188: LOGICAL LSAME
189: DOUBLE PRECISION DLAMCH, ZLANGE
190: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
191: EXTERNAL LSAME, DLAMCH, ZLANGE, ZDOTC, ZDOTU, ZLADIV
192: * ..
193: * .. External Subroutines ..
194: EXTERNAL DLABAD, XERBLA, ZDSCAL
195: * ..
196: * .. Intrinsic Functions ..
197: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
198: * ..
199: * .. Executable Statements ..
200: *
201: * Decode and Test input parameters
202: *
203: NOTRNA = LSAME( TRANA, 'N' )
204: NOTRNB = LSAME( TRANB, 'N' )
205: *
206: INFO = 0
207: IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'C' ) ) THEN
208: INFO = -1
209: ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'C' ) ) THEN
210: INFO = -2
211: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
212: INFO = -3
213: ELSE IF( M.LT.0 ) THEN
214: INFO = -4
215: ELSE IF( N.LT.0 ) THEN
216: INFO = -5
217: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
218: INFO = -7
219: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
220: INFO = -9
221: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
222: INFO = -11
223: END IF
224: IF( INFO.NE.0 ) THEN
225: CALL XERBLA( 'ZTRSYL', -INFO )
226: RETURN
227: END IF
228: *
229: * Quick return if possible
230: *
231: SCALE = ONE
232: IF( M.EQ.0 .OR. N.EQ.0 )
233: $ RETURN
234: *
235: * Set constants to control overflow
236: *
237: EPS = DLAMCH( 'P' )
238: SMLNUM = DLAMCH( 'S' )
239: BIGNUM = ONE / SMLNUM
240: CALL DLABAD( SMLNUM, BIGNUM )
241: SMLNUM = SMLNUM*DBLE( M*N ) / EPS
242: BIGNUM = ONE / SMLNUM
243: SMIN = MAX( SMLNUM, EPS*ZLANGE( 'M', M, M, A, LDA, DUM ),
244: $ EPS*ZLANGE( 'M', N, N, B, LDB, DUM ) )
245: SGN = ISGN
246: *
247: IF( NOTRNA .AND. NOTRNB ) THEN
248: *
249: * Solve A*X + ISGN*X*B = scale*C.
250: *
251: * The (K,L)th block of X is determined starting from
252: * bottom-left corner column by column by
253: *
254: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
255: *
256: * Where
257: * M L-1
258: * R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)].
259: * I=K+1 J=1
260: *
261: DO 30 L = 1, N
262: DO 20 K = M, 1, -1
263: *
264: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
265: $ C( MIN( K+1, M ), L ), 1 )
266: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
267: VEC = C( K, L ) - ( SUML+SGN*SUMR )
268: *
269: SCALOC = ONE
270: A11 = A( K, K ) + SGN*B( L, L )
271: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
272: IF( DA11.LE.SMIN ) THEN
273: A11 = SMIN
274: DA11 = SMIN
275: INFO = 1
276: END IF
277: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
278: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
279: IF( DB.GT.BIGNUM*DA11 )
280: $ SCALOC = ONE / DB
281: END IF
282: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
283: *
284: IF( SCALOC.NE.ONE ) THEN
285: DO 10 J = 1, N
286: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
287: 10 CONTINUE
288: SCALE = SCALE*SCALOC
289: END IF
290: C( K, L ) = X11
291: *
292: 20 CONTINUE
293: 30 CONTINUE
294: *
295: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
296: *
297: * Solve A**H *X + ISGN*X*B = scale*C.
298: *
299: * The (K,L)th block of X is determined starting from
300: * upper-left corner column by column by
301: *
302: * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
303: *
304: * Where
305: * K-1 L-1
306: * R(K,L) = SUM [A**H(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]
307: * I=1 J=1
308: *
309: DO 60 L = 1, N
310: DO 50 K = 1, M
311: *
312: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
313: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
314: VEC = C( K, L ) - ( SUML+SGN*SUMR )
315: *
316: SCALOC = ONE
317: A11 = DCONJG( A( K, K ) ) + SGN*B( L, L )
318: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
319: IF( DA11.LE.SMIN ) THEN
320: A11 = SMIN
321: DA11 = SMIN
322: INFO = 1
323: END IF
324: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
325: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
326: IF( DB.GT.BIGNUM*DA11 )
327: $ SCALOC = ONE / DB
328: END IF
329: *
330: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
331: *
332: IF( SCALOC.NE.ONE ) THEN
333: DO 40 J = 1, N
334: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
335: 40 CONTINUE
336: SCALE = SCALE*SCALOC
337: END IF
338: C( K, L ) = X11
339: *
340: 50 CONTINUE
341: 60 CONTINUE
342: *
343: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
344: *
345: * Solve A**H*X + ISGN*X*B**H = C.
346: *
347: * The (K,L)th block of X is determined starting from
348: * upper-right corner column by column by
349: *
350: * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
351: *
352: * Where
353: * K-1
354: * R(K,L) = SUM [A**H(I,K)*X(I,L)] +
355: * I=1
356: * N
357: * ISGN*SUM [X(K,J)*B**H(L,J)].
358: * J=L+1
359: *
360: DO 90 L = N, 1, -1
361: DO 80 K = 1, M
362: *
363: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
364: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
365: $ B( L, MIN( L+1, N ) ), LDB )
366: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
367: *
368: SCALOC = ONE
369: A11 = DCONJG( A( K, K )+SGN*B( L, L ) )
370: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
371: IF( DA11.LE.SMIN ) THEN
372: A11 = SMIN
373: DA11 = SMIN
374: INFO = 1
375: END IF
376: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
377: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
378: IF( DB.GT.BIGNUM*DA11 )
379: $ SCALOC = ONE / DB
380: END IF
381: *
382: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
383: *
384: IF( SCALOC.NE.ONE ) THEN
385: DO 70 J = 1, N
386: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
387: 70 CONTINUE
388: SCALE = SCALE*SCALOC
389: END IF
390: C( K, L ) = X11
391: *
392: 80 CONTINUE
393: 90 CONTINUE
394: *
395: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
396: *
397: * Solve A*X + ISGN*X*B**H = C.
398: *
399: * The (K,L)th block of X is determined starting from
400: * bottom-left corner column by column by
401: *
402: * A(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
403: *
404: * Where
405: * M N
406: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B**H(L,J)]
407: * I=K+1 J=L+1
408: *
409: DO 120 L = N, 1, -1
410: DO 110 K = M, 1, -1
411: *
412: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
413: $ C( MIN( K+1, M ), L ), 1 )
414: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
415: $ B( L, MIN( L+1, N ) ), LDB )
416: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
417: *
418: SCALOC = ONE
419: A11 = A( K, K ) + SGN*DCONJG( B( L, L ) )
420: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
421: IF( DA11.LE.SMIN ) THEN
422: A11 = SMIN
423: DA11 = SMIN
424: INFO = 1
425: END IF
426: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
427: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
428: IF( DB.GT.BIGNUM*DA11 )
429: $ SCALOC = ONE / DB
430: END IF
431: *
432: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
433: *
434: IF( SCALOC.NE.ONE ) THEN
435: DO 100 J = 1, N
436: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
437: 100 CONTINUE
438: SCALE = SCALE*SCALOC
439: END IF
440: C( K, L ) = X11
441: *
442: 110 CONTINUE
443: 120 CONTINUE
444: *
445: END IF
446: *
447: RETURN
448: *
449: * End of ZTRSYL
450: *
451: END
CVSweb interface <joel.bertrand@systella.fr>
![[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