294 SUBROUTINE ztgsyl( 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,
307 DOUBLE PRECISION DIF, SCALE
311 COMPLEX*16 A( lda, * ), B( ldb, * ), C( ldc, * ),
312 $ d( ldd, * ), e( lde, * ), f( ldf, * ),
321 DOUBLE PRECISION ZERO, ONE
322 parameter ( zero = 0.0d+0, one = 1.0d+0 )
324 parameter ( czero = (0.0d+0, 0.0d+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 DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC
335 EXTERNAL lsame, ilaenv
341 INTRINSIC dble, dcmplx, max, 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(
'ZTGSYL', -info )
398 ELSE IF( lquery )
THEN
404 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
416 mb = ilaenv( 2,
'ZTGSYL', trans, m, n, -1, -1 )
417 nb = ilaenv( 5,
'ZTGSYL', trans, m, n, -1, -1 )
424 CALL zlaset(
'F', m, n, czero, czero, c, ldc )
425 CALL zlaset(
'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 ztgsy2( 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( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
449 dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
452 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
457 CALL zlacpy(
'F', m, n, c, ldc, work, m )
458 CALL zlacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
459 CALL zlaset(
'F', m, n, czero, czero, c, ldc )
460 CALL zlaset(
'F', m, n, czero, czero, f, ldf )
461 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
462 CALL zlacpy(
'F', m, n, work, m, c, ldc )
463 CALL zlacpy(
'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 ztgsy2( 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 zscal( m, dcmplx( scaloc, zero ),
542 CALL zscal( m, dcmplx( scaloc, zero ),
546 CALL zscal( is-1, dcmplx( scaloc, zero ),
548 CALL zscal( is-1, dcmplx( scaloc, zero ),
552 CALL zscal( m-ie, dcmplx( scaloc, zero ),
554 CALL zscal( m-ie, dcmplx( scaloc, zero ),
558 CALL zscal( m, dcmplx( scaloc, zero ),
560 CALL zscal( m, dcmplx( scaloc, zero ),
569 CALL zgemm(
'N',
'N', is-1, nb, mb,
570 $ dcmplx( -one, zero ), a( 1, is ), lda,
571 $ c( is, js ), ldc, dcmplx( one, zero ),
573 CALL zgemm(
'N',
'N', is-1, nb, mb,
574 $ dcmplx( -one, zero ), d( 1, is ), ldd,
575 $ c( is, js ), ldc, dcmplx( one, zero ),
579 CALL zgemm(
'N',
'N', mb, n-je, nb,
580 $ dcmplx( one, zero ), f( is, js ), ldf,
581 $ b( js, je+1 ), ldb,
582 $ dcmplx( one, zero ), c( is, je+1 ),
584 CALL zgemm(
'N',
'N', mb, n-je, nb,
585 $ dcmplx( one, zero ), f( is, js ), ldf,
586 $ e( js, je+1 ), lde,
587 $ dcmplx( one, zero ), f( is, je+1 ),
592 IF( dscale.NE.zero )
THEN
593 IF( ijob.EQ.1 .OR. ijob.EQ.3 )
THEN
594 dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
596 dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
599 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
604 CALL zlacpy(
'F', m, n, c, ldc, work, m )
605 CALL zlacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
606 CALL zlaset(
'F', m, n, czero, czero, c, ldc )
607 CALL zlaset(
'F', m, n, czero, czero, f, ldf )
608 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
609 CALL zlacpy(
'F', m, n, work, m, c, ldc )
610 CALL zlacpy(
'F', m, n, work( m*n+1 ), m, f, ldf )
624 ie = iwork( i+1 ) - 1
626 DO 200 j = q, p + 2, -1
628 je = iwork( j+1 ) - 1
630 CALL ztgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
631 $ b( js, js ), ldb, c( is, js ), ldc,
632 $ d( is, is ), ldd, e( js, js ), lde,
633 $ f( is, js ), ldf, scaloc, dsum, dscale,
637 IF( scaloc.NE.one )
THEN
639 CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
641 CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
645 CALL zscal( is-1, dcmplx( scaloc, zero ),
647 CALL zscal( is-1, dcmplx( scaloc, zero ),
651 CALL zscal( m-ie, dcmplx( scaloc, zero ),
653 CALL zscal( m-ie, dcmplx( scaloc, zero ),
657 CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
659 CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
668 CALL zgemm(
'N',
'C', mb, js-1, nb,
669 $ dcmplx( one, zero ), c( is, js ), ldc,
670 $ b( 1, js ), ldb, dcmplx( one, zero ),
672 CALL zgemm(
'N',
'C', mb, js-1, nb,
673 $ dcmplx( one, zero ), f( is, js ), ldf,
674 $ e( 1, js ), lde, dcmplx( one, zero ),
678 CALL zgemm(
'C',
'N', m-ie, nb, mb,
679 $ dcmplx( -one, zero ), a( is, ie+1 ), lda,
680 $ c( is, js ), ldc, dcmplx( one, zero ),
681 $ c( ie+1, js ), ldc )
682 CALL zgemm(
'C',
'N', m-ie, nb, mb,
683 $ dcmplx( -one, zero ), d( is, ie+1 ), ldd,
684 $ f( is, js ), ldf, dcmplx( one, zero ),
685 $ c( ie+1, js ), ldc )
logical function lde(RI, RJ, LR)
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine ztgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, INFO)
ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine ztgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
ZTGSYL
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL