Annotation of rpl/lapack/lapack/ztrsyl.f, revision 1.9
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