001:       SUBROUTINE CTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
002:      $                   LDC, SCALE, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          TRANA, TRANB
011:       INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
012:       REAL               SCALE
013: *     ..
014: *     .. Array Arguments ..
015:       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CTRSYL solves the complex Sylvester matrix equation:
022: *
023: *     op(A)*X + X*op(B) = scale*C or
024: *     op(A)*X - X*op(B) = scale*C,
025: *
026: *  where op(A) = A or A**H, and A and B are both upper triangular. A is
027: *  M-by-M and B is N-by-N; the right hand side C and the solution X are
028: *  M-by-N; and scale is an output scale factor, set <= 1 to avoid
029: *  overflow in X.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  TRANA   (input) CHARACTER*1
035: *          Specifies the option op(A):
036: *          = 'N': op(A) = A    (No transpose)
037: *          = 'C': op(A) = A**H (Conjugate transpose)
038: *
039: *  TRANB   (input) CHARACTER*1
040: *          Specifies the option op(B):
041: *          = 'N': op(B) = B    (No transpose)
042: *          = 'C': op(B) = B**H (Conjugate transpose)
043: *
044: *  ISGN    (input) INTEGER
045: *          Specifies the sign in the equation:
046: *          = +1: solve op(A)*X + X*op(B) = scale*C
047: *          = -1: solve op(A)*X - X*op(B) = scale*C
048: *
049: *  M       (input) INTEGER
050: *          The order of the matrix A, and the number of rows in the
051: *          matrices X and C. M >= 0.
052: *
053: *  N       (input) INTEGER
054: *          The order of the matrix B, and the number of columns in the
055: *          matrices X and C. N >= 0.
056: *
057: *  A       (input) COMPLEX array, dimension (LDA,M)
058: *          The upper triangular matrix A.
059: *
060: *  LDA     (input) INTEGER
061: *          The leading dimension of the array A. LDA >= max(1,M).
062: *
063: *  B       (input) COMPLEX array, dimension (LDB,N)
064: *          The upper triangular matrix B.
065: *
066: *  LDB     (input) INTEGER
067: *          The leading dimension of the array B. LDB >= max(1,N).
068: *
069: *  C       (input/output) COMPLEX array, dimension (LDC,N)
070: *          On entry, the M-by-N right hand side matrix C.
071: *          On exit, C is overwritten by the solution matrix X.
072: *
073: *  LDC     (input) INTEGER
074: *          The leading dimension of the array C. LDC >= max(1,M)
075: *
076: *  SCALE   (output) REAL
077: *          The scale factor, scale, set <= 1 to avoid overflow in X.
078: *
079: *  INFO    (output) INTEGER
080: *          = 0: successful exit
081: *          < 0: if INFO = -i, the i-th argument had an illegal value
082: *          = 1: A and B have common or very close eigenvalues; perturbed
083: *               values were used to solve the equation (but the matrices
084: *               A and B are unchanged).
085: *
086: *  =====================================================================
087: *
088: *     .. Parameters ..
089:       REAL               ONE
090:       PARAMETER          ( ONE = 1.0E+0 )
091: *     ..
092: *     .. Local Scalars ..
093:       LOGICAL            NOTRNA, NOTRNB
094:       INTEGER            J, K, L
095:       REAL               BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
096:      $                   SMLNUM
097:       COMPLEX            A11, SUML, SUMR, VEC, X11
098: *     ..
099: *     .. Local Arrays ..
100:       REAL               DUM( 1 )
101: *     ..
102: *     .. External Functions ..
103:       LOGICAL            LSAME
104:       REAL               CLANGE, SLAMCH
105:       COMPLEX            CDOTC, CDOTU, CLADIV
106:       EXTERNAL           LSAME, CLANGE, SLAMCH, CDOTC, CDOTU, CLADIV
107: *     ..
108: *     .. External Subroutines ..
109:       EXTERNAL           CSSCAL, SLABAD, XERBLA
110: *     ..
111: *     .. Intrinsic Functions ..
112:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
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( 'CTRSYL', -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 = SLAMCH( 'P' )
153:       SMLNUM = SLAMCH( 'S' )
154:       BIGNUM = ONE / SMLNUM
155:       CALL SLABAD( SMLNUM, BIGNUM )
156:       SMLNUM = SMLNUM*REAL( M*N ) / EPS
157:       BIGNUM = ONE / SMLNUM
158:       SMIN = MAX( SMLNUM, EPS*CLANGE( 'M', M, M, A, LDA, DUM ),
159:      $       EPS*CLANGE( '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 = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
180:      $                C( MIN( K+1, M ), L ), 1 )
181:                SUMR = CDOTU( 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( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
187:                IF( DA11.LE.SMIN ) THEN
188:                   A11 = SMIN
189:                   DA11 = SMIN
190:                   INFO = 1
191:                END IF
192:                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( 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 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
198: *
199:                IF( SCALOC.NE.ONE ) THEN
200:                   DO 10 J = 1, N
201:                      CALL CSSCAL( 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 = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
228:                SUMR = CDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
229:                VEC = C( K, L ) - ( SUML+SGN*SUMR )
230: *
231:                SCALOC = ONE
232:                A11 = CONJG( A( K, K ) ) + SGN*B( L, L )
233:                DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
234:                IF( DA11.LE.SMIN ) THEN
235:                   A11 = SMIN
236:                   DA11 = SMIN
237:                   INFO = 1
238:                END IF
239:                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( 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 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
246: *
247:                IF( SCALOC.NE.ONE ) THEN
248:                   DO 40 J = 1, N
249:                      CALL CSSCAL( 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 = CDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
279:                SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
280:      $                B( L, MIN( L+1, N ) ), LDB )
281:                VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) )
282: *
283:                SCALOC = ONE
284:                A11 = CONJG( A( K, K )+SGN*B( L, L ) )
285:                DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
286:                IF( DA11.LE.SMIN ) THEN
287:                   A11 = SMIN
288:                   DA11 = SMIN
289:                   INFO = 1
290:                END IF
291:                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( 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 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
298: *
299:                IF( SCALOC.NE.ONE ) THEN
300:                   DO 70 J = 1, N
301:                      CALL CSSCAL( 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 = CDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
328:      $                C( MIN( K+1, M ), L ), 1 )
329:                SUMR = CDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
330:      $                B( L, MIN( L+1, N ) ), LDB )
331:                VEC = C( K, L ) - ( SUML+SGN*CONJG( SUMR ) )
332: *
333:                SCALOC = ONE
334:                A11 = A( K, K ) + SGN*CONJG( B( L, L ) )
335:                DA11 = ABS( REAL( A11 ) ) + ABS( AIMAG( A11 ) )
336:                IF( DA11.LE.SMIN ) THEN
337:                   A11 = SMIN
338:                   DA11 = SMIN
339:                   INFO = 1
340:                END IF
341:                DB = ABS( REAL( VEC ) ) + ABS( AIMAG( 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 = CLADIV( VEC*CMPLX( SCALOC ), A11 )
348: *
349:                IF( SCALOC.NE.ONE ) THEN
350:                   DO 100 J = 1, N
351:                      CALL CSSCAL( 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 CTRSYL
365: *
366:       END
367: