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, * )
300 parameter ( ldz = 8 )
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 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
logical function lde(RI, RJ, LR)
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sgetc2(N, A, LDA, IPIV, JPIV, INFO)
SGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
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 sgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
SGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
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 saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY