294 SUBROUTINE ctgsyl( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
295 $ ldd, e,
lde, f, ldf, scale, dif, work, lwork,
305 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
311 COMPLEX A( lda, * ), B( ldb, * ), C( ldc, * ),
312 $ d( ldd, * ), e( lde, * ), f( ldf, * ),
322 parameter ( zero = 0.0e+0, one = 1.0e+0 )
324 parameter ( czero = (0.0e+0, 0.0e+0) )
327 LOGICAL LQUERY, NOTRAN
328 INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
329 $ linfo, lwmin, mb, nb, p, pq, q
330 REAL DSCALE, DSUM, SCALE2, SCALOC
335 EXTERNAL lsame, ilaenv
341 INTRINSIC cmplx, max,
REAL, SQRT
348 notran = lsame( trans,
'N' )
349 lquery = ( lwork.EQ.-1 )
351 IF( .NOT.notran .AND. .NOT.lsame( trans,
'C' ) )
THEN
353 ELSE IF( notran )
THEN
354 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) )
THEN
361 ELSE IF( n.LE.0 )
THEN
363 ELSE IF( lda.LT.max( 1, m ) )
THEN
365 ELSE IF( ldb.LT.max( 1, n ) )
THEN
367 ELSE IF( ldc.LT.max( 1, m ) )
THEN
369 ELSE IF( ldd.LT.max( 1, m ) )
THEN
371 ELSE IF( lde.LT.max( 1, n ) )
THEN
373 ELSE IF( ldf.LT.max( 1, m ) )
THEN
380 IF( ijob.EQ.1 .OR. ijob.EQ.2 )
THEN
381 lwmin = max( 1, 2*m*n )
390 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
396 CALL xerbla(
'CTGSYL', -info )
398 ELSE IF( lquery )
THEN
404 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
416 mb = ilaenv( 2,
'CTGSYL', trans, m, n, -1, -1 )
417 nb = ilaenv( 5,
'CTGSYL', trans, m, n, -1, -1 )
424 CALL claset(
'F', m, n, czero, czero, c, ldc )
425 CALL claset(
'F', m, n, czero, czero, f, ldf )
426 ELSE IF( ijob.GE.1 .AND. notran )
THEN
431 IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
436 DO 30 iround = 1, isolve
442 CALL ctgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
443 $ ldd, e, lde, f, ldf, scale, dsum, dscale,
445 IF( dscale.NE.zero )
THEN
446 IF( ijob.EQ.1 .OR. ijob.EQ.3 )
THEN
447 dif = sqrt(
REAL( 2*M*N ) ) / ( DSCALE*SQRT( dsum ) )
449 dif = sqrt(
REAL( PQ ) ) / ( DSCALE*SQRT( dsum ) )
452 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
457 CALL clacpy(
'F', m, n, c, ldc, work, m )
458 CALL clacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
459 CALL claset(
'F', m, n, czero, czero, c, ldc )
460 CALL claset(
'F', m, n, czero, czero, f, ldf )
461 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
462 CALL clacpy(
'F', m, n, work, m, c, ldc )
463 CALL clacpy(
'F', m, n, work( m*n+1 ), m, f, ldf )
487 IF( iwork( p ).EQ.iwork( p+1 ) )
507 IF( iwork( q ).EQ.iwork( q+1 ) )
511 DO 150 iround = 1, isolve
524 je = iwork( j+1 ) - 1
528 ie = iwork( i+1 ) - 1
530 CALL ctgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
531 $ b( js, js ), ldb, c( is, js ), ldc,
532 $ d( is, is ), ldd, e( js, js ), lde,
533 $ f( is, js ), ldf, scaloc, dsum, dscale,
538 IF( scaloc.NE.one )
THEN
540 CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
542 CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
546 CALL cscal( is-1, cmplx( scaloc, zero ),
548 CALL cscal( is-1, cmplx( scaloc, zero ),
552 CALL cscal( m-ie, cmplx( scaloc, zero ),
554 CALL cscal( m-ie, cmplx( scaloc, zero ),
558 CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
560 CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
569 CALL cgemm(
'N',
'N', is-1, nb, mb,
570 $ cmplx( -one, zero ), a( 1, is ), lda,
571 $ c( is, js ), ldc, cmplx( one, zero ),
573 CALL cgemm(
'N',
'N', is-1, nb, mb,
574 $ cmplx( -one, zero ), d( 1, is ), ldd,
575 $ c( is, js ), ldc, cmplx( one, zero ),
579 CALL cgemm(
'N',
'N', mb, n-je, nb,
580 $ cmplx( one, zero ), f( is, js ), ldf,
581 $ b( js, je+1 ), ldb, cmplx( one, zero ),
582 $ c( is, je+1 ), ldc )
583 CALL cgemm(
'N',
'N', mb, n-je, nb,
584 $ cmplx( one, zero ), f( is, js ), ldf,
585 $ e( js, je+1 ), lde, cmplx( one, zero ),
586 $ f( is, je+1 ), ldf )
590 IF( dscale.NE.zero )
THEN
591 IF( ijob.EQ.1 .OR. ijob.EQ.3 )
THEN
592 dif = sqrt(
REAL( 2*M*N ) ) / ( DSCALE*SQRT( dsum ) )
594 dif = sqrt(
REAL( PQ ) ) / ( DSCALE*SQRT( dsum ) )
597 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
602 CALL clacpy(
'F', m, n, c, ldc, work, m )
603 CALL clacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
604 CALL claset(
'F', m, n, czero, czero, c, ldc )
605 CALL claset(
'F', m, n, czero, czero, f, ldf )
606 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
607 CALL clacpy(
'F', m, n, work, m, c, ldc )
608 CALL clacpy(
'F', m, n, work( m*n+1 ), m, f, ldf )
622 ie = iwork( i+1 ) - 1
624 DO 200 j = q, p + 2, -1
626 je = iwork( j+1 ) - 1
628 CALL ctgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
629 $ b( js, js ), ldb, c( is, js ), ldc,
630 $ d( is, is ), ldd, e( js, js ), lde,
631 $ f( is, js ), ldf, scaloc, dsum, dscale,
635 IF( scaloc.NE.one )
THEN
637 CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
639 CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
643 CALL cscal( is-1, cmplx( scaloc, zero ), c( 1, k ),
645 CALL cscal( is-1, cmplx( scaloc, zero ), f( 1, k ),
649 CALL cscal( m-ie, cmplx( scaloc, zero ),
651 CALL cscal( m-ie, cmplx( scaloc, zero ),
655 CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
657 CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
666 CALL cgemm(
'N',
'C', mb, js-1, nb,
667 $ cmplx( one, zero ), c( is, js ), ldc,
668 $ b( 1, js ), ldb, cmplx( one, zero ),
670 CALL cgemm(
'N',
'C', mb, js-1, nb,
671 $ cmplx( one, zero ), f( is, js ), ldf,
672 $ e( 1, js ), lde, cmplx( one, zero ),
676 CALL cgemm(
'C',
'N', m-ie, nb, mb,
677 $ cmplx( -one, zero ), a( is, ie+1 ), lda,
678 $ c( is, js ), ldc, cmplx( one, zero ),
679 $ c( ie+1, js ), ldc )
680 CALL cgemm(
'C',
'N', m-ie, nb, mb,
681 $ cmplx( -one, zero ), d( is, ie+1 ), ldd,
682 $ f( is, js ), ldf, cmplx( one, zero ),
683 $ c( ie+1, js ), ldc )
logical function lde(RI, RJ, LR)
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine ctgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, INFO)
CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
subroutine ctgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
CTGSYL
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM