273 SUBROUTINE stgsy2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
274 $ ldd, e, lde, f, ldf, scale, rdsum, rdscal,
284 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
286 REAL RDSCAL, RDSUM, SCALE
290 REAL A( lda, * ), B( ldb, * ), C( ldc, * ),
291 $ d( ldd, * ), e( lde, * ), f( ldf, * )
302 parameter( zero = 0.0e+0, one = 1.0e+0 )
306 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
307 $ k, mb, nb, p, q, zdim
311 INTEGER IPIV( ldz ), JPIV( ldz )
312 REAL RHS( ldz ), Z( ldz, ldz )
331 notran = lsame( trans,
'N' )
332 IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) )
THEN
334 ELSE IF( notran )
THEN
335 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) )
THEN
342 ELSE IF( n.LE.0 )
THEN
344 ELSE IF( lda.LT.max( 1, m ) )
THEN
346 ELSE IF( ldb.LT.max( 1, n ) )
THEN
348 ELSE IF( ldc.LT.max( 1, m ) )
THEN
350 ELSE IF( ldd.LT.max( 1, m ) )
THEN
352 ELSE IF( lde.LT.max( 1, n ) )
THEN
354 ELSE IF( ldf.LT.max( 1, m ) )
THEN
359 CALL
xerbla(
'STGSY2', -info )
375 IF( a( i+1, i ).NE.zero )
THEN
395 IF( b( j+1, j ).NE.zero )
THEN
417 je = iwork( j+1 ) - 1
423 ie = iwork( i+1 ) - 1
427 IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) )
THEN
431 z( 1, 1 ) = a( is, is )
432 z( 2, 1 ) = d( is, is )
433 z( 1, 2 ) = -b( js, js )
434 z( 2, 2 ) = -e( js, js )
438 rhs( 1 ) = c( is, js )
439 rhs( 2 ) = f( is, js )
443 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
448 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
450 IF( scaloc.NE.one )
THEN
452 CALL
sscal( m, scaloc, c( 1, k ), 1 )
453 CALL
sscal( m, scaloc, f( 1, k ), 1 )
458 CALL
slatdf( ijob, zdim, z, ldz, rhs, rdsum,
459 $ rdscal, ipiv, jpiv )
464 c( is, js ) = rhs( 1 )
465 f( is, js ) = rhs( 2 )
472 CALL
saxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
474 CALL
saxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
478 CALL
saxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
479 $ c( is, je+1 ), ldc )
480 CALL
saxpy( n-je, rhs( 2 ), e( js, je+1 ), lde,
481 $ f( is, je+1 ), ldf )
484 ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) )
THEN
488 z( 1, 1 ) = a( is, is )
490 z( 3, 1 ) = d( is, is )
494 z( 2, 2 ) = a( is, is )
496 z( 4, 2 ) = d( is, is )
498 z( 1, 3 ) = -b( js, js )
499 z( 2, 3 ) = -b( js, jsp1 )
500 z( 3, 3 ) = -e( js, js )
501 z( 4, 3 ) = -e( js, jsp1 )
503 z( 1, 4 ) = -b( jsp1, js )
504 z( 2, 4 ) = -b( jsp1, jsp1 )
506 z( 4, 4 ) = -e( jsp1, jsp1 )
510 rhs( 1 ) = c( is, js )
511 rhs( 2 ) = c( is, jsp1 )
512 rhs( 3 ) = f( is, js )
513 rhs( 4 ) = f( is, jsp1 )
517 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
522 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
524 IF( scaloc.NE.one )
THEN
526 CALL
sscal( m, scaloc, c( 1, k ), 1 )
527 CALL
sscal( m, scaloc, f( 1, k ), 1 )
532 CALL
slatdf( ijob, zdim, z, ldz, rhs, rdsum,
533 $ rdscal, ipiv, jpiv )
538 c( is, js ) = rhs( 1 )
539 c( is, jsp1 ) = rhs( 2 )
540 f( is, js ) = rhs( 3 )
541 f( is, jsp1 ) = rhs( 4 )
547 CALL
sger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
548 $ 1, c( 1, js ), ldc )
549 CALL
sger( is-1, nb, -one, d( 1, is ), 1, rhs( 1 ),
550 $ 1, f( 1, js ), ldf )
553 CALL
saxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
554 $ c( is, je+1 ), ldc )
555 CALL
saxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
556 $ f( is, je+1 ), ldf )
557 CALL
saxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
558 $ c( is, je+1 ), ldc )
559 CALL
saxpy( n-je, rhs( 4 ), e( jsp1, je+1 ), lde,
560 $ f( is, je+1 ), ldf )
563 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) )
THEN
567 z( 1, 1 ) = a( is, is )
568 z( 2, 1 ) = a( isp1, is )
569 z( 3, 1 ) = d( is, is )
572 z( 1, 2 ) = a( is, isp1 )
573 z( 2, 2 ) = a( isp1, isp1 )
574 z( 3, 2 ) = d( is, isp1 )
575 z( 4, 2 ) = d( isp1, isp1 )
577 z( 1, 3 ) = -b( js, js )
579 z( 3, 3 ) = -e( js, js )
583 z( 2, 4 ) = -b( js, js )
585 z( 4, 4 ) = -e( js, js )
589 rhs( 1 ) = c( is, js )
590 rhs( 2 ) = c( isp1, js )
591 rhs( 3 ) = f( is, js )
592 rhs( 4 ) = f( isp1, js )
596 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
600 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
602 IF( scaloc.NE.one )
THEN
604 CALL
sscal( m, scaloc, c( 1, k ), 1 )
605 CALL
sscal( m, scaloc, f( 1, k ), 1 )
610 CALL
slatdf( ijob, zdim, z, ldz, rhs, rdsum,
611 $ rdscal, ipiv, jpiv )
616 c( is, js ) = rhs( 1 )
617 c( isp1, js ) = rhs( 2 )
618 f( is, js ) = rhs( 3 )
619 f( isp1, js ) = rhs( 4 )
625 CALL
sgemv(
'N', is-1, mb, -one, a( 1, is ), lda,
626 $ rhs( 1 ), 1, one, c( 1, js ), 1 )
627 CALL
sgemv(
'N', is-1, mb, -one, d( 1, is ), ldd,
628 $ rhs( 1 ), 1, one, f( 1, js ), 1 )
631 CALL
sger( mb, n-je, one, rhs( 3 ), 1,
632 $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
633 CALL
sger( mb, n-je, one, rhs( 3 ), 1,
634 $ e( js, je+1 ), lde, f( is, je+1 ), ldf )
637 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) )
THEN
641 CALL
slaset(
'F', ldz, ldz, zero, zero, z, ldz )
643 z( 1, 1 ) = a( is, is )
644 z( 2, 1 ) = a( isp1, is )
645 z( 5, 1 ) = d( is, is )
647 z( 1, 2 ) = a( is, isp1 )
648 z( 2, 2 ) = a( isp1, isp1 )
649 z( 5, 2 ) = d( is, isp1 )
650 z( 6, 2 ) = d( isp1, isp1 )
652 z( 3, 3 ) = a( is, is )
653 z( 4, 3 ) = a( isp1, is )
654 z( 7, 3 ) = d( is, is )
656 z( 3, 4 ) = a( is, isp1 )
657 z( 4, 4 ) = a( isp1, isp1 )
658 z( 7, 4 ) = d( is, isp1 )
659 z( 8, 4 ) = d( isp1, isp1 )
661 z( 1, 5 ) = -b( js, js )
662 z( 3, 5 ) = -b( js, jsp1 )
663 z( 5, 5 ) = -e( js, js )
664 z( 7, 5 ) = -e( js, jsp1 )
666 z( 2, 6 ) = -b( js, js )
667 z( 4, 6 ) = -b( js, jsp1 )
668 z( 6, 6 ) = -e( js, js )
669 z( 8, 6 ) = -e( js, jsp1 )
671 z( 1, 7 ) = -b( jsp1, js )
672 z( 3, 7 ) = -b( jsp1, jsp1 )
673 z( 7, 7 ) = -e( jsp1, jsp1 )
675 z( 2, 8 ) = -b( jsp1, js )
676 z( 4, 8 ) = -b( jsp1, jsp1 )
677 z( 8, 8 ) = -e( jsp1, jsp1 )
684 CALL
scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
685 CALL
scopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
692 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
696 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
698 IF( scaloc.NE.one )
THEN
700 CALL
sscal( m, scaloc, c( 1, k ), 1 )
701 CALL
sscal( m, scaloc, f( 1, k ), 1 )
706 CALL
slatdf( ijob, zdim, z, ldz, rhs, rdsum,
707 $ rdscal, ipiv, jpiv )
714 DO 100 jj = 0, nb - 1
715 CALL
scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
716 CALL
scopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
725 CALL
sgemm(
'N',
'N', is-1, nb, mb, -one,
726 $ a( 1, is ), lda, rhs( 1 ), mb, one,
728 CALL
sgemm(
'N',
'N', is-1, nb, mb, -one,
729 $ d( 1, is ), ldd, rhs( 1 ), mb, one,
734 CALL
sgemm(
'N',
'N', mb, n-je, nb, one, rhs( k ),
735 $ mb, b( js, je+1 ), ldb, one,
736 $ c( is, je+1 ), ldc )
737 CALL
sgemm(
'N',
'N', mb, n-je, nb, one, rhs( k ),
738 $ mb, e( js, je+1 ), lde, one,
739 $ f( is, je+1 ), ldf )
759 ie = iwork( i+1 ) - 1
761 DO 190 j = q, p + 2, -1
765 je = iwork( j+1 ) - 1
768 IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) )
THEN
772 z( 1, 1 ) = a( is, is )
773 z( 2, 1 ) = -b( js, js )
774 z( 1, 2 ) = d( is, is )
775 z( 2, 2 ) = -e( js, js )
779 rhs( 1 ) = c( is, js )
780 rhs( 2 ) = f( is, js )
784 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
788 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
789 IF( scaloc.NE.one )
THEN
791 CALL
sscal( m, scaloc, c( 1, k ), 1 )
792 CALL
sscal( m, scaloc, f( 1, k ), 1 )
799 c( is, js ) = rhs( 1 )
800 f( is, js ) = rhs( 2 )
807 CALL
saxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
810 CALL
saxpy( js-1, alpha, e( 1, js ), 1, f( is, 1 ),
815 CALL
saxpy( m-ie, alpha, a( is, ie+1 ), lda,
818 CALL
saxpy( m-ie, alpha, d( is, ie+1 ), ldd,
822 ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) )
THEN
826 z( 1, 1 ) = a( is, is )
828 z( 3, 1 ) = -b( js, js )
829 z( 4, 1 ) = -b( jsp1, js )
832 z( 2, 2 ) = a( is, is )
833 z( 3, 2 ) = -b( js, jsp1 )
834 z( 4, 2 ) = -b( jsp1, jsp1 )
836 z( 1, 3 ) = d( is, is )
838 z( 3, 3 ) = -e( js, js )
842 z( 2, 4 ) = d( is, is )
843 z( 3, 4 ) = -e( js, jsp1 )
844 z( 4, 4 ) = -e( jsp1, jsp1 )
848 rhs( 1 ) = c( is, js )
849 rhs( 2 ) = c( is, jsp1 )
850 rhs( 3 ) = f( is, js )
851 rhs( 4 ) = f( is, jsp1 )
855 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
858 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
859 IF( scaloc.NE.one )
THEN
861 CALL
sscal( m, scaloc, c( 1, k ), 1 )
862 CALL
sscal( m, scaloc, f( 1, k ), 1 )
869 c( is, js ) = rhs( 1 )
870 c( is, jsp1 ) = rhs( 2 )
871 f( is, js ) = rhs( 3 )
872 f( is, jsp1 ) = rhs( 4 )
878 CALL
saxpy( js-1, rhs( 1 ), b( 1, js ), 1,
880 CALL
saxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
882 CALL
saxpy( js-1, rhs( 3 ), e( 1, js ), 1,
884 CALL
saxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
888 CALL
sger( m-ie, nb, -one, a( is, ie+1 ), lda,
889 $ rhs( 1 ), 1, c( ie+1, js ), ldc )
890 CALL
sger( m-ie, nb, -one, d( is, ie+1 ), ldd,
891 $ rhs( 3 ), 1, c( ie+1, js ), ldc )
894 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) )
THEN
898 z( 1, 1 ) = a( is, is )
899 z( 2, 1 ) = a( is, isp1 )
900 z( 3, 1 ) = -b( js, js )
903 z( 1, 2 ) = a( isp1, is )
904 z( 2, 2 ) = a( isp1, isp1 )
906 z( 4, 2 ) = -b( js, js )
908 z( 1, 3 ) = d( is, is )
909 z( 2, 3 ) = d( is, isp1 )
910 z( 3, 3 ) = -e( js, js )
914 z( 2, 4 ) = d( isp1, isp1 )
916 z( 4, 4 ) = -e( js, js )
920 rhs( 1 ) = c( is, js )
921 rhs( 2 ) = c( isp1, js )
922 rhs( 3 ) = f( is, js )
923 rhs( 4 ) = f( isp1, js )
927 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
931 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
932 IF( scaloc.NE.one )
THEN
934 CALL
sscal( m, scaloc, c( 1, k ), 1 )
935 CALL
sscal( m, scaloc, f( 1, k ), 1 )
942 c( is, js ) = rhs( 1 )
943 c( isp1, js ) = rhs( 2 )
944 f( is, js ) = rhs( 3 )
945 f( isp1, js ) = rhs( 4 )
951 CALL
sger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
952 $ 1, f( is, 1 ), ldf )
953 CALL
sger( mb, js-1, one, rhs( 3 ), 1, e( 1, js ),
954 $ 1, f( is, 1 ), ldf )
957 CALL
sgemv(
'T', mb, m-ie, -one, a( is, ie+1 ),
958 $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
960 CALL
sgemv(
'T', mb, m-ie, -one, d( is, ie+1 ),
961 $ ldd, rhs( 3 ), 1, one, c( ie+1, js ),
965 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) )
THEN
969 CALL
slaset(
'F', ldz, ldz, zero, zero, z, ldz )
971 z( 1, 1 ) = a( is, is )
972 z( 2, 1 ) = a( is, isp1 )
973 z( 5, 1 ) = -b( js, js )
974 z( 7, 1 ) = -b( jsp1, js )
976 z( 1, 2 ) = a( isp1, is )
977 z( 2, 2 ) = a( isp1, isp1 )
978 z( 6, 2 ) = -b( js, js )
979 z( 8, 2 ) = -b( jsp1, js )
981 z( 3, 3 ) = a( is, is )
982 z( 4, 3 ) = a( is, isp1 )
983 z( 5, 3 ) = -b( js, jsp1 )
984 z( 7, 3 ) = -b( jsp1, jsp1 )
986 z( 3, 4 ) = a( isp1, is )
987 z( 4, 4 ) = a( isp1, isp1 )
988 z( 6, 4 ) = -b( js, jsp1 )
989 z( 8, 4 ) = -b( jsp1, jsp1 )
991 z( 1, 5 ) = d( is, is )
992 z( 2, 5 ) = d( is, isp1 )
993 z( 5, 5 ) = -e( js, js )
995 z( 2, 6 ) = d( isp1, isp1 )
996 z( 6, 6 ) = -e( js, js )
998 z( 3, 7 ) = d( is, is )
999 z( 4, 7 ) = d( is, isp1 )
1000 z( 5, 7 ) = -e( js, jsp1 )
1001 z( 7, 7 ) = -e( jsp1, jsp1 )
1003 z( 4, 8 ) = d( isp1, isp1 )
1004 z( 6, 8 ) = -e( js, jsp1 )
1005 z( 8, 8 ) = -e( jsp1, jsp1 )
1011 DO 160 jj = 0, nb - 1
1012 CALL
scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1013 CALL
scopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
1021 CALL
sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1025 CALL
sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1026 IF( scaloc.NE.one )
THEN
1028 CALL
sscal( m, scaloc, c( 1, k ), 1 )
1029 CALL
sscal( m, scaloc, f( 1, k ), 1 )
1031 scale = scale*scaloc
1038 DO 180 jj = 0, nb - 1
1039 CALL
scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1040 CALL
scopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
1049 CALL
sgemm(
'N',
'T', mb, js-1, nb, one,
1050 $ c( is, js ), ldc, b( 1, js ), ldb, one,
1052 CALL
sgemm(
'N',
'T', mb, js-1, nb, one,
1053 $ f( is, js ), ldf, e( 1, js ), lde, one,
1057 CALL
sgemm(
'T',
'N', m-ie, nb, mb, -one,
1058 $ a( is, ie+1 ), lda, c( is, js ), ldc,
1059 $ one, c( ie+1, js ), ldc )
1060 CALL
sgemm(
'T',
'N', m-ie, nb, mb, -one,
1061 $ d( is, ie+1 ), ldd, f( is, js ), ldf,
1062 $ one, c( ie+1, js ), ldc )
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine stgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
STGSY2 solves the generalized Sylvester equation (unblocked algorithm).
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sgetc2(N, A, LDA, IPIV, JPIV, INFO)
SGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
SGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
subroutine sscal(N, SA, SX, INCX)
SSCAL