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, * )
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 )
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 dgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
LOGICAL function lsame(CA, CB)
LSAME
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
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 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 ...
LOGICAL function lde(RI, RJ, LR)
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dgetc2(N, A, LDA, IPIV, JPIV, INFO)
DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...
subroutine xerbla(SRNAME, INFO)
XERBLA
for consistency with libblas dev *Do not ship a copy of jquery js
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dscal(N, DA, DX, INCX)
DSCAL
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
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 dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV