273 SUBROUTINE dtgsy2( 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 DOUBLE PRECISION RDSCAL, RDSUM, SCALE
290 DOUBLE PRECISION A( lda, * ), B( ldb, * ), C( ldc, * ),
291 $ d( ldd, * ), e( lde, * ), f( ldf, * )
300 parameter ( ldz = 8 )
301 DOUBLE PRECISION ZERO, ONE
302 parameter ( zero = 0.0d+0, one = 1.0d+0 )
306 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
307 $ k, mb, nb, p, q, zdim
308 DOUBLE PRECISION ALPHA, SCALOC
311 INTEGER IPIV( ldz ), JPIV( ldz )
312 DOUBLE PRECISION 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(
'DTGSY2', -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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
448 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
450 IF( scaloc.NE.one )
THEN
452 CALL dscal( m, scaloc, c( 1, k ), 1 )
453 CALL dscal( m, scaloc, f( 1, k ), 1 )
458 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
459 $ rdscal, ipiv, jpiv )
464 c( is, js ) = rhs( 1 )
465 f( is, js ) = rhs( 2 )
472 CALL daxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
474 CALL daxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
478 CALL daxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
479 $ c( is, je+1 ), ldc )
480 CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
522 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
524 IF( scaloc.NE.one )
THEN
526 CALL dscal( m, scaloc, c( 1, k ), 1 )
527 CALL dscal( m, scaloc, f( 1, k ), 1 )
532 CALL dlatdf( 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 dger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
548 $ 1, c( 1, js ), ldc )
549 CALL dger( is-1, nb, -one, d( 1, is ), 1, rhs( 1 ),
550 $ 1, f( 1, js ), ldf )
553 CALL daxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
554 $ c( is, je+1 ), ldc )
555 CALL daxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
556 $ f( is, je+1 ), ldf )
557 CALL daxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
558 $ c( is, je+1 ), ldc )
559 CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
600 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
602 IF( scaloc.NE.one )
THEN
604 CALL dscal( m, scaloc, c( 1, k ), 1 )
605 CALL dscal( m, scaloc, f( 1, k ), 1 )
610 CALL dlatdf( 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 dgemv(
'N', is-1, mb, -one, a( 1, is ), lda,
626 $ rhs( 1 ), 1, one, c( 1, js ), 1 )
627 CALL dgemv(
'N', is-1, mb, -one, d( 1, is ), ldd,
628 $ rhs( 1 ), 1, one, f( 1, js ), 1 )
631 CALL dger( mb, n-je, one, rhs( 3 ), 1,
632 $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
633 CALL dger( 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 dlaset(
'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 dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
685 CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
692 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
696 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
698 IF( scaloc.NE.one )
THEN
700 CALL dscal( m, scaloc, c( 1, k ), 1 )
701 CALL dscal( m, scaloc, f( 1, k ), 1 )
706 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
707 $ rdscal, ipiv, jpiv )
714 DO 100 jj = 0, nb - 1
715 CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
716 CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
725 CALL dgemm(
'N',
'N', is-1, nb, mb, -one,
726 $ a( 1, is ), lda, rhs( 1 ), mb, one,
728 CALL dgemm(
'N',
'N', is-1, nb, mb, -one,
729 $ d( 1, is ), ldd, rhs( 1 ), mb, one,
734 CALL dgemm(
'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 dgemm(
'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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
788 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
789 IF( scaloc.NE.one )
THEN
791 CALL dscal( m, scaloc, c( 1, k ), 1 )
792 CALL dscal( m, scaloc, f( 1, k ), 1 )
799 c( is, js ) = rhs( 1 )
800 f( is, js ) = rhs( 2 )
807 CALL daxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
810 CALL daxpy( js-1, alpha, e( 1, js ), 1, f( is, 1 ),
815 CALL daxpy( m-ie, alpha, a( is, ie+1 ), lda,
818 CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
858 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
859 IF( scaloc.NE.one )
THEN
861 CALL dscal( m, scaloc, c( 1, k ), 1 )
862 CALL dscal( 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 daxpy( js-1, rhs( 1 ), b( 1, js ), 1,
880 CALL daxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
882 CALL daxpy( js-1, rhs( 3 ), e( 1, js ), 1,
884 CALL daxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
888 CALL dger( m-ie, nb, -one, a( is, ie+1 ), lda,
889 $ rhs( 1 ), 1, c( ie+1, js ), ldc )
890 CALL dger( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
931 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
932 IF( scaloc.NE.one )
THEN
934 CALL dscal( m, scaloc, c( 1, k ), 1 )
935 CALL dscal( 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 dger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
952 $ 1, f( is, 1 ), ldf )
953 CALL dger( mb, js-1, one, rhs( 3 ), 1, e( 1, js ),
954 $ 1, f( is, 1 ), ldf )
957 CALL dgemv(
'T', mb, m-ie, -one, a( is, ie+1 ),
958 $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
960 CALL dgemv(
'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 dlaset(
'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 dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1013 CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
1021 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1025 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1026 IF( scaloc.NE.one )
THEN
1028 CALL dscal( m, scaloc, c( 1, k ), 1 )
1029 CALL dscal( m, scaloc, f( 1, k ), 1 )
1031 scale = scale*scaloc
1038 DO 180 jj = 0, nb - 1
1039 CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1040 CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
1049 CALL dgemm(
'N',
'T', mb, js-1, nb, one,
1050 $ c( is, js ), ldc, b( 1, js ), ldb, one,
1052 CALL dgemm(
'N',
'T', mb, js-1, nb, one,
1053 $ f( is, js ), ldf, e( 1, js ), lde, one,
1057 CALL dgemm(
'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 dgemm(
'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 dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
logical function lde(RI, RJ, LR)
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dtgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
subroutine dgetc2(N, A, LDA, IPIV, JPIV, INFO)
DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...