LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctrsyl ( character  TRANA,
character  TRANB,
integer  ISGN,
integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldc, * )  C,
integer  LDC,
real  SCALE,
integer  INFO 
)

CTRSYL

Download CTRSYL + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 CTRSYL solves the complex Sylvester matrix equation:

    op(A)*X + X*op(B) = scale*C or
    op(A)*X - X*op(B) = scale*C,

 where op(A) = A or A**H, and A and B are both upper triangular. A is
 M-by-M and B is N-by-N; the right hand side C and the solution X are
 M-by-N; and scale is an output scale factor, set <= 1 to avoid
 overflow in X.
Parameters
[in]TRANA
          TRANA is CHARACTER*1
          Specifies the option op(A):
          = 'N': op(A) = A    (No transpose)
          = 'C': op(A) = A**H (Conjugate transpose)
[in]TRANB
          TRANB is CHARACTER*1
          Specifies the option op(B):
          = 'N': op(B) = B    (No transpose)
          = 'C': op(B) = B**H (Conjugate transpose)
[in]ISGN
          ISGN is INTEGER
          Specifies the sign in the equation:
          = +1: solve op(A)*X + X*op(B) = scale*C
          = -1: solve op(A)*X - X*op(B) = scale*C
[in]M
          M is INTEGER
          The order of the matrix A, and the number of rows in the
          matrices X and C. M >= 0.
[in]N
          N is INTEGER
          The order of the matrix B, and the number of columns in the
          matrices X and C. N >= 0.
[in]A
          A is COMPLEX array, dimension (LDA,M)
          The upper triangular matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in]B
          B is COMPLEX array, dimension (LDB,N)
          The upper triangular matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[in,out]C
          C is COMPLEX array, dimension (LDC,N)
          On entry, the M-by-N right hand side matrix C.
          On exit, C is overwritten by the solution matrix X.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M)
[out]SCALE
          SCALE is REAL
          The scale factor, scale, set <= 1 to avoid overflow in X.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          = 1: A and B have common or very close eigenvalues; perturbed
               values were used to solve the equation (but the matrices
               A and B are unchanged).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 159 of file ctrsyl.f.

159 *
160 * -- LAPACK computational routine (version 3.4.0) --
161 * -- LAPACK is a software package provided by Univ. of Tennessee, --
162 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163 * November 2011
164 *
165 * .. Scalar Arguments ..
166  CHARACTER trana, tranb
167  INTEGER info, isgn, lda, ldb, ldc, m, n
168  REAL scale
169 * ..
170 * .. Array Arguments ..
171  COMPLEX a( lda, * ), b( ldb, * ), c( ldc, * )
172 * ..
173 *
174 * =====================================================================
175 *
176 * .. Parameters ..
177  REAL one
178  parameter ( one = 1.0e+0 )
179 * ..
180 * .. Local Scalars ..
181  LOGICAL notrna, notrnb
182  INTEGER j, k, l
183  REAL bignum, da11, db, eps, scaloc, sgn, smin,
184  $ smlnum
185  COMPLEX a11, suml, sumr, vec, x11
186 * ..
187 * .. Local Arrays ..
188  REAL dum( 1 )
189 * ..
190 * .. External Functions ..
191  LOGICAL lsame
192  REAL clange, slamch
193  COMPLEX cdotc, cdotu, cladiv
194  EXTERNAL lsame, clange, slamch, cdotc, cdotu, cladiv
195 * ..
196 * .. External Subroutines ..
197  EXTERNAL csscal, slabad, xerbla
198 * ..
199 * .. Intrinsic Functions ..
200  INTRINSIC abs, aimag, cmplx, conjg, max, min, real
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( 'CTRSYL', -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 = slamch( 'P' )
241  smlnum = slamch( 'S' )
242  bignum = one / smlnum
243  CALL slabad( smlnum, bignum )
244  smlnum = smlnum*REAL( M*N ) / eps
245  bignum = one / smlnum
246  smin = max( smlnum, eps*clange( 'M', m, m, a, lda, dum ),
247  $ eps*clange( '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 = cdotu( m-k, a( k, min( k+1, m ) ), lda,
268  $ c( min( k+1, m ), l ), 1 )
269  sumr = cdotu( 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( REAL( A11 ) ) + abs( aimag( a11 ) )
275  IF( da11.LE.smin ) THEN
276  a11 = smin
277  da11 = smin
278  info = 1
279  END IF
280  db = abs( REAL( VEC ) ) + abs( aimag( 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 = cladiv( vec*cmplx( scaloc ), a11 )
286 *
287  IF( scaloc.NE.one ) THEN
288  DO 10 j = 1, n
289  CALL csscal( 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 *
300 * Solve A**H *X + ISGN*X*B = scale*C.
301 *
302 * The (K,L)th block of X is determined starting from
303 * upper-left corner column by column by
304 *
305 * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
306 *
307 * Where
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
311 *
312  DO 60 l = 1, n
313  DO 50 k = 1, m
314 *
315  suml = cdotc( k-1, a( 1, k ), 1, c( 1, l ), 1 )
316  sumr = cdotu( l-1, c( k, 1 ), ldc, b( 1, l ), 1 )
317  vec = c( k, l ) - ( suml+sgn*sumr )
318 *
319  scaloc = one
320  a11 = conjg( a( k, k ) ) + sgn*b( l, l )
321  da11 = abs( REAL( A11 ) ) + abs( aimag( a11 ) )
322  IF( da11.LE.smin ) THEN
323  a11 = smin
324  da11 = smin
325  info = 1
326  END IF
327  db = abs( REAL( VEC ) ) + abs( aimag( 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 = cladiv( vec*cmplx( scaloc ), a11 )
334 *
335  IF( scaloc.NE.one ) THEN
336  DO 40 j = 1, n
337  CALL csscal( 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 *
348 * Solve A**H*X + ISGN*X*B**H = C.
349 *
350 * The (K,L)th block of X is determined starting from
351 * upper-right corner column by column by
352 *
353 * A**H(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
354 *
355 * Where
356 * K-1
357 * R(K,L) = SUM [A**H(I,K)*X(I,L)] +
358 * I=1
359 * N
360 * ISGN*SUM [X(K,J)*B**H(L,J)].
361 * J=L+1
362 *
363  DO 90 l = n, 1, -1
364  DO 80 k = 1, m
365 *
366  suml = cdotc( k-1, a( 1, k ), 1, c( 1, l ), 1 )
367  sumr = cdotc( n-l, c( k, min( l+1, n ) ), ldc,
368  $ b( l, min( l+1, n ) ), ldb )
369  vec = c( k, l ) - ( suml+sgn*conjg( sumr ) )
370 *
371  scaloc = one
372  a11 = conjg( a( k, k )+sgn*b( l, l ) )
373  da11 = abs( REAL( A11 ) ) + abs( aimag( a11 ) )
374  IF( da11.LE.smin ) THEN
375  a11 = smin
376  da11 = smin
377  info = 1
378  END IF
379  db = abs( REAL( VEC ) ) + abs( aimag( 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 = cladiv( vec*cmplx( scaloc ), a11 )
386 *
387  IF( scaloc.NE.one ) THEN
388  DO 70 j = 1, n
389  CALL csscal( 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 *
400 * Solve A*X + ISGN*X*B**H = C.
401 *
402 * The (K,L)th block of X is determined starting from
403 * bottom-left corner column by column by
404 *
405 * A(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
406 *
407 * Where
408 * M N
409 * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B**H(L,J)]
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 = cdotu( m-k, a( k, min( k+1, m ) ), lda,
416  $ c( min( k+1, m ), l ), 1 )
417  sumr = cdotc( n-l, c( k, min( l+1, n ) ), ldc,
418  $ b( l, min( l+1, n ) ), ldb )
419  vec = c( k, l ) - ( suml+sgn*conjg( sumr ) )
420 *
421  scaloc = one
422  a11 = a( k, k ) + sgn*conjg( b( l, l ) )
423  da11 = abs( REAL( A11 ) ) + abs( aimag( a11 ) )
424  IF( da11.LE.smin ) THEN
425  a11 = smin
426  da11 = smin
427  info = 1
428  END IF
429  db = abs( REAL( VEC ) ) + abs( aimag( 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 = cladiv( vec*cmplx( scaloc ), a11 )
436 *
437  IF( scaloc.NE.one ) THEN
438  DO 100 j = 1, n
439  CALL csscal( 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 CTRSYL
453 *
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
clange
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cdotu
complex function cdotu(N, CX, INCX, CY, INCY)
CDOTU
Definition: cdotu.f:54
cdotc
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:54
cladiv
complex function cladiv(X, Y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: cladiv.f:66
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: