Annotation of rpl/lapack/lapack/ztrsyl.f, revision 1.14
1.9 bertrand 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: *> \date November 2011
153: *
154: *> \ingroup complex16SYcomputational
155: *
156: * =====================================================================
1.1 bertrand 157: SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
158: $ LDC, SCALE, INFO )
159: *
1.9 bertrand 160: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 161: * -- LAPACK is a software package provided by Univ. of Tennessee, --
162: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 bertrand 163: * November 2011
1.1 bertrand 164: *
165: * .. Scalar Arguments ..
166: CHARACTER TRANA, TRANB
167: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
168: DOUBLE PRECISION SCALE
169: * ..
170: * .. Array Arguments ..
171: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
172: * ..
173: *
174: * =====================================================================
175: *
176: * .. Parameters ..
177: DOUBLE PRECISION ONE
178: PARAMETER ( ONE = 1.0D+0 )
179: * ..
180: * .. Local Scalars ..
181: LOGICAL NOTRNA, NOTRNB
182: INTEGER J, K, L
183: DOUBLE PRECISION BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
184: $ SMLNUM
185: COMPLEX*16 A11, SUML, SUMR, VEC, X11
186: * ..
187: * .. Local Arrays ..
188: DOUBLE PRECISION DUM( 1 )
189: * ..
190: * .. External Functions ..
191: LOGICAL LSAME
192: DOUBLE PRECISION DLAMCH, ZLANGE
193: COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
194: EXTERNAL LSAME, DLAMCH, ZLANGE, ZDOTC, ZDOTU, ZLADIV
195: * ..
196: * .. External Subroutines ..
197: EXTERNAL DLABAD, XERBLA, ZDSCAL
198: * ..
199: * .. Intrinsic Functions ..
200: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
201: * ..
202: * .. Executable Statements ..
203: *
204: * Decode and Test input parameters
205: *
206: NOTRNA = LSAME( TRANA, 'N' )
207: NOTRNB = LSAME( TRANB, 'N' )
208: *
209: INFO = 0
210: IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'C' ) ) THEN
211: INFO = -1
212: ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'C' ) ) THEN
213: INFO = -2
214: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
215: INFO = -3
216: ELSE IF( M.LT.0 ) THEN
217: INFO = -4
218: ELSE IF( N.LT.0 ) THEN
219: INFO = -5
220: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
221: INFO = -7
222: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
223: INFO = -9
224: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
225: INFO = -11
226: END IF
227: IF( INFO.NE.0 ) THEN
228: CALL XERBLA( 'ZTRSYL', -INFO )
229: RETURN
230: END IF
231: *
232: * Quick return if possible
233: *
234: SCALE = ONE
235: IF( M.EQ.0 .OR. N.EQ.0 )
236: $ RETURN
237: *
238: * Set constants to control overflow
239: *
240: EPS = DLAMCH( 'P' )
241: SMLNUM = DLAMCH( 'S' )
242: BIGNUM = ONE / SMLNUM
243: CALL DLABAD( SMLNUM, BIGNUM )
244: SMLNUM = SMLNUM*DBLE( M*N ) / EPS
245: BIGNUM = ONE / SMLNUM
246: SMIN = MAX( SMLNUM, EPS*ZLANGE( 'M', M, M, A, LDA, DUM ),
247: $ EPS*ZLANGE( 'M', N, N, B, LDB, DUM ) )
248: SGN = ISGN
249: *
250: IF( NOTRNA .AND. NOTRNB ) THEN
251: *
252: * Solve A*X + ISGN*X*B = scale*C.
253: *
254: * The (K,L)th block of X is determined starting from
255: * bottom-left corner column by column by
256: *
257: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
258: *
259: * Where
260: * M L-1
261: * R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)].
262: * I=K+1 J=1
263: *
264: DO 30 L = 1, N
265: DO 20 K = M, 1, -1
266: *
267: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
268: $ C( MIN( K+1, M ), L ), 1 )
269: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
270: VEC = C( K, L ) - ( SUML+SGN*SUMR )
271: *
272: SCALOC = ONE
273: A11 = A( K, K ) + SGN*B( L, L )
274: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
275: IF( DA11.LE.SMIN ) THEN
276: A11 = SMIN
277: DA11 = SMIN
278: INFO = 1
279: END IF
280: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
281: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
282: IF( DB.GT.BIGNUM*DA11 )
283: $ SCALOC = ONE / DB
284: END IF
285: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
286: *
287: IF( SCALOC.NE.ONE ) THEN
288: DO 10 J = 1, N
289: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
290: 10 CONTINUE
291: SCALE = SCALE*SCALOC
292: END IF
293: C( K, L ) = X11
294: *
295: 20 CONTINUE
296: 30 CONTINUE
297: *
298: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
299: *
1.8 bertrand 300: * Solve A**H *X + ISGN*X*B = scale*C.
1.1 bertrand 301: *
302: * The (K,L)th block of X is determined starting from
303: * upper-left corner column by column by
304: *
1.8 bertrand 305: * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
1.1 bertrand 306: *
307: * Where
1.8 bertrand 308: * K-1 L-1
309: * R(K,L) = SUM [A**H(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]
310: * I=1 J=1
1.1 bertrand 311: *
312: DO 60 L = 1, N
313: DO 50 K = 1, M
314: *
315: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
316: SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
317: VEC = C( K, L ) - ( SUML+SGN*SUMR )
318: *
319: SCALOC = ONE
320: A11 = DCONJG( A( K, K ) ) + SGN*B( L, L )
321: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
322: IF( DA11.LE.SMIN ) THEN
323: A11 = SMIN
324: DA11 = SMIN
325: INFO = 1
326: END IF
327: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
328: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
329: IF( DB.GT.BIGNUM*DA11 )
330: $ SCALOC = ONE / DB
331: END IF
332: *
333: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
334: *
335: IF( SCALOC.NE.ONE ) THEN
336: DO 40 J = 1, N
337: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
338: 40 CONTINUE
339: SCALE = SCALE*SCALOC
340: END IF
341: C( K, L ) = X11
342: *
343: 50 CONTINUE
344: 60 CONTINUE
345: *
346: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
347: *
1.8 bertrand 348: * Solve A**H*X + ISGN*X*B**H = C.
1.1 bertrand 349: *
350: * The (K,L)th block of X is determined starting from
351: * upper-right corner column by column by
352: *
1.8 bertrand 353: * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
1.1 bertrand 354: *
355: * Where
356: * K-1
1.8 bertrand 357: * R(K,L) = SUM [A**H(I,K)*X(I,L)] +
1.1 bertrand 358: * I=1
359: * N
1.8 bertrand 360: * ISGN*SUM [X(K,J)*B**H(L,J)].
1.1 bertrand 361: * J=L+1
362: *
363: DO 90 L = N, 1, -1
364: DO 80 K = 1, M
365: *
366: SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
367: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
368: $ B( L, MIN( L+1, N ) ), LDB )
369: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
370: *
371: SCALOC = ONE
372: A11 = DCONJG( A( K, K )+SGN*B( L, L ) )
373: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
374: IF( DA11.LE.SMIN ) THEN
375: A11 = SMIN
376: DA11 = SMIN
377: INFO = 1
378: END IF
379: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
380: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
381: IF( DB.GT.BIGNUM*DA11 )
382: $ SCALOC = ONE / DB
383: END IF
384: *
385: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
386: *
387: IF( SCALOC.NE.ONE ) THEN
388: DO 70 J = 1, N
389: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
390: 70 CONTINUE
391: SCALE = SCALE*SCALOC
392: END IF
393: C( K, L ) = X11
394: *
395: 80 CONTINUE
396: 90 CONTINUE
397: *
398: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
399: *
1.8 bertrand 400: * Solve A*X + ISGN*X*B**H = C.
1.1 bertrand 401: *
402: * The (K,L)th block of X is determined starting from
403: * bottom-left corner column by column by
404: *
1.8 bertrand 405: * A(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
1.1 bertrand 406: *
407: * Where
408: * M N
1.8 bertrand 409: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B**H(L,J)]
1.1 bertrand 410: * I=K+1 J=L+1
411: *
412: DO 120 L = N, 1, -1
413: DO 110 K = M, 1, -1
414: *
415: SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
416: $ C( MIN( K+1, M ), L ), 1 )
417: SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
418: $ B( L, MIN( L+1, N ) ), LDB )
419: VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
420: *
421: SCALOC = ONE
422: A11 = A( K, K ) + SGN*DCONJG( B( L, L ) )
423: DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
424: IF( DA11.LE.SMIN ) THEN
425: A11 = SMIN
426: DA11 = SMIN
427: INFO = 1
428: END IF
429: DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
430: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
431: IF( DB.GT.BIGNUM*DA11 )
432: $ SCALOC = ONE / DB
433: END IF
434: *
435: X11 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
436: *
437: IF( SCALOC.NE.ONE ) THEN
438: DO 100 J = 1, N
439: CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
440: 100 CONTINUE
441: SCALE = SCALE*SCALOC
442: END IF
443: C( K, L ) = X11
444: *
445: 110 CONTINUE
446: 120 CONTINUE
447: *
448: END IF
449: *
450: RETURN
451: *
452: * End of ZTRSYL
453: *
454: 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