1       SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
  2      $                   LDC, SCALE, INFO )
  3 *
  4 *  -- LAPACK routine (version 3.3.1) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *  -- April 2011                                                      --
  8 *
  9 *     .. Scalar Arguments ..
 10       CHARACTER          TRANA, TRANB
 11       INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
 12       DOUBLE PRECISION   SCALE
 13 *     ..
 14 *     .. Array Arguments ..
 15       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  ZTRSYL solves the complex Sylvester matrix equation:
 22 *
 23 *     op(A)*X + X*op(B) = scale*C or
 24 *     op(A)*X - X*op(B) = scale*C,
 25 *
 26 *  where op(A) = A or A**H, and A and B are both upper triangular. A is
 27 *  M-by-M and B is N-by-N; the right hand side C and the solution X are
 28 *  M-by-N; and scale is an output scale factor, set <= 1 to avoid
 29 *  overflow in X.
 30 *
 31 *  Arguments
 32 *  =========
 33 *
 34 *  TRANA   (input) CHARACTER*1
 35 *          Specifies the option op(A):
 36 *          = 'N': op(A) = A    (No transpose)
 37 *          = 'C': op(A) = A**H (Conjugate transpose)
 38 *
 39 *  TRANB   (input) CHARACTER*1
 40 *          Specifies the option op(B):
 41 *          = 'N': op(B) = B    (No transpose)
 42 *          = 'C': op(B) = B**H (Conjugate transpose)
 43 *
 44 *  ISGN    (input) INTEGER
 45 *          Specifies the sign in the equation:
 46 *          = +1: solve op(A)*X + X*op(B) = scale*C
 47 *          = -1: solve op(A)*X - X*op(B) = scale*C
 48 *
 49 *  M       (input) INTEGER
 50 *          The order of the matrix A, and the number of rows in the
 51 *          matrices X and C. M >= 0.
 52 *
 53 *  N       (input) INTEGER
 54 *          The order of the matrix B, and the number of columns in the
 55 *          matrices X and C. N >= 0.
 56 *
 57 *  A       (input) COMPLEX*16 array, dimension (LDA,M)
 58 *          The upper triangular matrix A.
 59 *
 60 *  LDA     (input) INTEGER
 61 *          The leading dimension of the array A. LDA >= max(1,M).
 62 *
 63 *  B       (input) COMPLEX*16 array, dimension (LDB,N)
 64 *          The upper triangular matrix B.
 65 *
 66 *  LDB     (input) INTEGER
 67 *          The leading dimension of the array B. LDB >= max(1,N).
 68 *
 69 *  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
 70 *          On entry, the M-by-N right hand side matrix C.
 71 *          On exit, C is overwritten by the solution matrix X.
 72 *
 73 *  LDC     (input) INTEGER
 74 *          The leading dimension of the array C. LDC >= max(1,M)
 75 *
 76 *  SCALE   (output) DOUBLE PRECISION
 77 *          The scale factor, scale, set <= 1 to avoid overflow in X.
 78 *
 79 *  INFO    (output) INTEGER
 80 *          = 0: successful exit
 81 *          < 0: if INFO = -i, the i-th argument had an illegal value
 82 *          = 1: A and B have common or very close eigenvalues; perturbed
 83 *               values were used to solve the equation (but the matrices
 84 *               A and B are unchanged).
 85 *
 86 *  =====================================================================
 87 *
 88 *     .. Parameters ..
 89       DOUBLE PRECISION   ONE
 90       PARAMETER          ( ONE = 1.0D+0 )
 91 *     ..
 92 *     .. Local Scalars ..
 93       LOGICAL            NOTRNA, NOTRNB
 94       INTEGER            J, K, L
 95       DOUBLE PRECISION   BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
 96      $                   SMLNUM
 97       COMPLEX*16         A11, SUML, SUMR, VEC, X11
 98 *     ..
 99 *     .. Local Arrays ..
100       DOUBLE PRECISION   DUM( 1 )
101 *     ..
102 *     .. External Functions ..
103       LOGICAL            LSAME
104       DOUBLE PRECISION   DLAMCH, ZLANGE
105       COMPLEX*16         ZDOTC, ZDOTU, ZLADIV
106       EXTERNAL           LSAME, DLAMCH, ZLANGE, ZDOTC, ZDOTU, ZLADIV
107 *     ..
108 *     .. External Subroutines ..
109       EXTERNAL           DLABAD, XERBLA, ZDSCAL
110 *     ..
111 *     .. Intrinsic Functions ..
112       INTRINSIC          ABSDBLEDCMPLXDCONJGDIMAGMAXMIN
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.MAX1, M ) ) THEN
133          INFO = -7
134       ELSE IF( LDB.LT.MAX1, N ) ) THEN
135          INFO = -9
136       ELSE IF( LDC.LT.MAX1, M ) ) THEN
137          INFO = -11
138       END IF
139       IF( INFO.NE.0 ) THEN
140          CALL XERBLA( 'ZTRSYL'-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 = DLAMCH( 'P' )
153       SMLNUM = DLAMCH( 'S' )
154       BIGNUM = ONE / SMLNUM
155       CALL DLABAD( SMLNUM, BIGNUM )
156       SMLNUM = SMLNUM*DBLE( M*N ) / EPS
157       BIGNUM = ONE / SMLNUM
158       SMIN = MAX( SMLNUM, EPS*ZLANGE( 'M', M, M, A, LDA, DUM ),
159      $       EPS*ZLANGE( '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 = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
180      $                C( MIN( K+1, M ), L ), 1 )
181                SUMR = ZDOTU( 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 = ABSDBLE( A11 ) ) + ABSDIMAG( A11 ) )
187                IF( DA11.LE.SMIN ) THEN
188                   A11 = SMIN
189                   DA11 = SMIN
190                   INFO = 1
191                END IF
192                DB = ABSDBLE( VEC ) ) + ABSDIMAG( 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 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
198 *
199                IF( SCALOC.NE.ONE ) THEN
200                   DO 10 J = 1, N
201                      CALL ZDSCAL( 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**H *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**H(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**H(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 = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
228                SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
229                VEC = C( K, L ) - ( SUML+SGN*SUMR )
230 *
231                SCALOC = ONE
232                A11 = DCONJG( A( K, K ) ) + SGN*B( L, L )
233                DA11 = ABSDBLE( A11 ) ) + ABSDIMAG( A11 ) )
234                IF( DA11.LE.SMIN ) THEN
235                   A11 = SMIN
236                   DA11 = SMIN
237                   INFO = 1
238                END IF
239                DB = ABSDBLE( VEC ) ) + ABSDIMAG( 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 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
246 *
247                IF( SCALOC.NE.ONE ) THEN
248                   DO 40 J = 1, N
249                      CALL ZDSCAL( 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**H*X + ISGN*X*B**H = 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**H(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
266 *
267 *        Where
268 *                    K-1
269 *           R(K,L) = SUM [A**H(I,K)*X(I,L)] +
270 *                    I=1
271 *                           N
272 *                     ISGN*SUM [X(K,J)*B**H(L,J)].
273 *                          J=L+1
274 *
275          DO 90 L = N, 1-1
276             DO 80 K = 1, M
277 *
278                SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
279                SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
280      $                B( L, MIN( L+1, N ) ), LDB )
281                VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
282 *
283                SCALOC = ONE
284                A11 = DCONJG( A( K, K )+SGN*B( L, L ) )
285                DA11 = ABSDBLE( A11 ) ) + ABSDIMAG( A11 ) )
286                IF( DA11.LE.SMIN ) THEN
287                   A11 = SMIN
288                   DA11 = SMIN
289                   INFO = 1
290                END IF
291                DB = ABSDBLE( VEC ) ) + ABSDIMAG( 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 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
298 *
299                IF( SCALOC.NE.ONE ) THEN
300                   DO 70 J = 1, N
301                      CALL ZDSCAL( 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**H = 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**H(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**H(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 = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
328      $                C( MIN( K+1, M ), L ), 1 )
329                SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
330      $                B( L, MIN( L+1, N ) ), LDB )
331                VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
332 *
333                SCALOC = ONE
334                A11 = A( K, K ) + SGN*DCONJG( B( L, L ) )
335                DA11 = ABSDBLE( A11 ) ) + ABSDIMAG( A11 ) )
336                IF( DA11.LE.SMIN ) THEN
337                   A11 = SMIN
338                   DA11 = SMIN
339                   INFO = 1
340                END IF
341                DB = ABSDBLE( VEC ) ) + ABSDIMAG( 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 = ZLADIV( VEC*DCMPLX( SCALOC ), A11 )
348 *
349                IF( SCALOC.NE.ONE ) THEN
350                   DO 100 J = 1, N
351                      CALL ZDSCAL( 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 ZTRSYL
365 *
366       END