3 USE v3_utilities,
ONLY:
assert
9 subroutine pdtrdf(m,nblock,Amat,Bmat,Cmat,ipivmat, desc)
29 REAL(rprec),
parameter :: one = 1.0d0
30 REAL(rprec),
parameter :: zero = 0.0d0
32 integer,
intent(in) :: m, nblock
33 REAL(rprec),
target,
dimension(:,:),
intent(inout) :: Amat
34 REAL(rprec),
target,
dimension(:,:),
intent(inout) :: Bmat
35 REAL(rprec),
target,
dimension(:,:),
intent(inout) :: Cmat
36 integer,
dimension(:,:),
target,
intent(inout) :: ipivmat
37 integer,
dimension(:),
intent(in) :: desc
40 integer :: k, ia,ja, ib,jb, iu,ju, ic,jc
41 integer :: info, nrhs, mm,nn
42 integer,
dimension(:),
pointer :: ipiv
43 REAL(rprec) :: alpha, beta
44 REAL(rprec),
dimension(:),
pointer :: Ak, Bkm1, Ukm1
45 REAL(rprec),
dimension(:),
pointer :: Lk, Ck
72 call pdgemm(
'N',
'N', m,m,m,
73 & alpha, bkm1,ib,jb,desc,
75 & beta, ak, ia,ja, desc )
88 call pdgetrf(mm,nn,lk,ia,ja,desc,ipiv,info)
89 call assert(info.eq.0,
'pdtrdf: pdgetrf return info != 0')
91 if (k .le. (nblock-1))
then
96 call pdgetrs(
'No',m,nrhs,lk,ia,ja,desc,ipiv,
97 & ck,ic,jc,desc, info)
98 call assert(info.eq.0,
'pdtrdf: pdgetrs return info != 0')
104 end subroutine pdtrdf
106 subroutine pdtrds(m,nblock,Amat,Bmat,Cmat,ipivmat,desc, &
107 & nrhs, Rrhs_in,ir,jr,descRrhs_in)
108 use descriptor_mod, ir1=>ir, jr1=>jr, mm1=>mm, info1=>info
123 integer,
parameter :: idebug = 0
124 REAL(rprec),
parameter :: one = 1.0d0
125 REAL(rprec),
parameter :: zero = 0.0d0
127 integer,
intent(in) :: m, nblock
128 REAL(rprec),
dimension(:,:),
target,
intent(inout) :: Amat
129 REAL(rprec),
dimension(:,:),
target,
intent(inout) :: Bmat
130 REAL(rprec),
dimension(:,:),
target,
intent(inout) :: Cmat
131 integer,
dimension(:,:),
target,
intent(inout) :: ipivmat
132 integer,
dimension(:),
intent(in) :: desc
134 integer,
intent(in) :: nrhs
135 integer,
intent(in) :: ir,jr
136 REAL(rprec),
dimension(:) :: Rrhs_in
137 integer,
dimension(:),
intent(in) :: descRrhs_in
140 integer :: k, info, mm,nn,kk
141 integer :: ia,ja, irk, jrk, iy,jy, ix,jx, iu,ju, ib,jb
142 REAL(rprec) :: alpha, beta
144 REAL(rprec),
pointer,
dimension(:) :: Lk,Bkm1,Uk
145 REAL(rprec),
pointer,
dimension(:) :: Rrhsk,Rrhskm1,Rrhskp1
146 integer,
pointer,
dimension(:) :: ipiv
149 integer :: inc1,inc2,iblock,irhs
151 logical,
parameter :: use_pdcopy = .false.
152 REAL(rprec),
dimension(:,:),
target,
allocatable :: Rrhs
153 integer,
dimension(DLEN_) :: descRrhs
161 icontxt = descrrhs_in(ctxt_)
162 mb = descrrhs_in(mb_)
167 call blacs_gridinfo( icontxt, nprow,npcol,myrow,mycol)
168 locp = numroc( m, mb, myrow,rsrc,nprow)
169 locq = numroc( nrhs,nb,mycol,csrc,npcol)
171 call descinit(descrrhs,m,nrhs,mb,nb,rsrc,csrc,icontxt,lld,info)
172 call assert(info.eq.0,
'pdtrds: descinit return info != 0')
174 ineed = max(1,locp*locq)
175 allocate( rrhs(ineed,nblock),stat=ierr)
176 call assert(ierr.eq.0,
'pdtrds: alloc Rrhs,ierr != 0')
191 ia = (ir-1) + 1 + (iblock-1)*m
199 rrhsk => rrhs(:,iblock)
200 call pdcopy(m, rrhs_in,ia,ja,descrrhs_in,inc1, &
201 & rrhsk,ib,jb,descrrhs,inc2 )
209 ia = (ir-1) + 1 + (iblock-1)*m
214 rrhsk => rrhs(:,iblock)
217 call pdgeadd(
'N',m,nrhs,alpha,rrhs_in,ia,ja,descrrhs_in, &
218 & beta, rrhsk,ib,jb,descrrhs )
233 icontxt = desc(ctxt_)
234 call blacs_gridinfo(icontxt, nprow,npcol,myrow,mycol)
235 isroot = (myrow.eq.0).and.(mycol.eq.0)
253 if (isroot .and. (idebug.ge.1))
then
254 write(*,*)
'pdtrds 77: m,nblock,nrhs ',m,nblock,nrhs
255 write(*,*)
'descRrhs_in(M_) ',descrrhs_in(m_)
256 write(*,*)
'descRrhs_in(N_) ',descrrhs_in(n_)
257 write(*,*)
'descRrhs_in(MB_) ',descrrhs_in(mb_)
258 write(*,*)
'descRrhs_in(NB_) ',descrrhs_in(nb_)
281 rrhskm1 => rrhs(:,k-1)
283 call pdgemm(
'N',
'N', mm,nn,kk,
284 & alpha, bkm1, ib,jb, desc,
285 & rrhskm1, iy,jy, descrrhs,
286 & beta, rrhsk, irk,jrk,descrrhs )
303 if (isroot .and. (idebug.ge.1))
then
304 write(*,*)
'pdtrds 106: k,irk,jrk ',k,irk,jrk
306 call pdgetrs(
'N', m,nrhs, lk,ia,ja,desc, ipiv,
307 & rrhsk,irk,jrk,descrrhs,info)
308 call assert(info.eq.0,
'pdtrds: pdgetrs return info != 0')
348 rrhskp1 => rrhs(:,k+1)
351 call pdgemm(
'N',
'N', mm,nn,kk,
352 & alpha, uk,iu,ju,desc,
353 & rrhskp1,ix,jx,descrrhs,
354 & beta, rrhsk,irk,jrk,descrrhs )
365 ia = (ir-1) + 1 + (iblock-1)*m
373 rrhsk => rrhs(:,iblock)
374 call pdcopy(m,rrhsk,ib,jb,descrrhs,inc2,
375 & rrhs_in,ia,ja,descrrhs_in,inc1)
383 ia = (ir-1) + 1 + (iblock-1)*m
388 rrhsk => rrhs(:,iblock)
392 call pdgeadd(
'N',m,nrhs,alpha,rrhsk,ib,jb,descrrhs, &
393 & beta,rrhs_in,ia,ja,descrrhs_in)
397 deallocate( rrhs, stat=ierr)
398 call assert(ierr.eq.0,
'pdtrds: dealloc Rhs ')
402 end subroutine pdtrds
404 #if defined(NEED_TOOLS)
405 SUBROUTINE pdgetrs( TRANS, N, NRHS, A, IA, JA, DESCA, IPIV, B, &
406 & IB, JB, DESCB, INFO )
415 INTEGER IA, IB, INFO, JA, JB, N, NRHS
418 INTEGER DESCA( * ), DESCB( * ), IPIV( * )
419 REAL(dp) A( * ), B( * )
556 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
559 & ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
560 & rsrc_ = 7, csrc_ = 8, lld_ = 9 )
562 PARAMETER ( ONE = 1.0d+0 )
566 INTEGER IAROW, IBROW, ICOFFA, ICTXT, IROFFA, IROFFB,
570 INTEGER DESCIP( DLEN_ ), IDUM1( 1 ), IDUM2( 1 )
573 EXTERNAL blacs_gridinfo, chk1mat, descset, pchk2mat,
574 & pdlapiv, pdtrsm, pxerbla
578 INTEGER INDXG2P, NUMROC
579 EXTERNAL indxg2p, lsame, numroc
588 ictxt = desca( ctxt_ )
589 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
594 IF( nprow.EQ.-1 )
THEN
597 notran = lsame( trans,
'N' )
598 CALL chk1mat( n, 2, n, 2, ia, ja, desca, 7, info )
599 CALL chk1mat( n, 2, nrhs, 3, ib, jb, descb, 12, info )
601 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
603 ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
605 iroffa = mod( ia-1, desca( mb_ ) )
606 icoffa = mod( ja-1, desca( nb_ ) )
607 iroffb = mod( ib-1, descb( mb_ ) )
608 IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) .AND. .NOT.
609 & lsame( trans,
'C' ) )
THEN
611 ELSE IF( iroffa.NE.0 )
THEN
613 ELSE IF( icoffa.NE.0 )
THEN
615 ELSE IF( desca( mb_ ).NE.desca( nb_ ) )
THEN
619 ELSE IF( descb( mb_ ).NE.desca( nb_ ) )
THEN
621 ELSE IF( ictxt.NE.descb( ctxt_ ) )
THEN
626 idum1( 1 ) = ichar(
'N' )
627 ELSE IF( lsame( trans,
'T' ) )
THEN
628 idum1( 1 ) = ichar(
'T' )
630 idum1( 1 ) = ichar(
'C' )
633 CALL pchk2mat( n, 2, n, 2, ia, ja, desca, 7, n, 2, nrhs, 3,
634 & ib, jb, descb, 12, 1, idum1, idum2, info )
638 CALL pxerbla( ictxt,
'PDGETRS', -info )
644 IF( n.EQ.0 .OR. nrhs.EQ.0 )
RETURN
646 CALL descset( descip, desca( m_ ) + desca( mb_ )*nprow, 1,
647 & desca( mb_ ), 1, desca( rsrc_ ), mycol, ictxt,
648 & desca( mb_ ) + numroc( desca( m_ ), desca( mb_ ),
649 & myrow, desca( rsrc_ ), nprow ) )
657 CALL pdlapiv(
'Forward',
'Row',
'Col', n, nrhs, b, ib, jb,
658 & descb, ipiv, ia, 1, descip, idum1 )
662 CALL pdtrsm(
'Left',
'Lower',
'No transpose',
'Unit', n, nrhs,
663 & one, a, ia, ja, desca, b, ib, jb, descb )
667 CALL pdtrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', n,
668 & nrhs, one, a, ia, ja, desca, b, ib, jb, descb )
675 CALL pdtrsm(
'Left',
'Upper',
'Transpose',
'Non-unit', n, nrhs,
676 & one, a, ia, ja, desca, b, ib, jb, descb )
680 CALL pdtrsm(
'Left',
'Lower',
'Transpose',
'Unit', n, nrhs,
681 & one, a, ia, ja, desca, b, ib, jb, descb )
685 CALL pdlapiv(
'Backward',
'Row',
'Col', n, nrhs, b, ib, jb,
686 & descb, ipiv, ia, 1, descip, idum1 )
694 END SUBROUTINE pdgetrs
696 SUBROUTINE pdelget( SCOPE, TOP, ALPHA, A, IA, JA, DESCA )
697 use descriptor_mod, desca_x=>desca
706 CHARACTER*1 SCOPE, TOP
818 PARAMETER ( ZERO = 0.0d+0 )
821 INTEGER IACOL, IAROW, ICTXT, IIA, IOFFA, JJA
824 EXTERNAL blacs_gridinfo, dgebr2d, dgebs2d, infog2l
834 ictxt = desca( ctxt_ )
835 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
837 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
842 IF( lsame( scope,
'R' ) )
THEN
843 IF( myrow.EQ.iarow )
THEN
844 IF( mycol.EQ.iacol )
THEN
845 ioffa = iia+(jja-1)*desca( lld_ )
846 CALL dgebs2d( ictxt, scope, top, 1, 1, a( ioffa ), 1 )
849 CALL dgebr2d( ictxt, scope, top, 1, 1, alpha, 1,
853 ELSE IF( lsame( scope,
'C' ) )
THEN
854 IF( mycol.EQ.iacol )
THEN
855 IF( myrow.EQ.iarow )
THEN
856 ioffa = iia+(jja-1)*desca( lld_ )
857 CALL dgebs2d( ictxt, scope, top, 1, 1, a( ioffa ), 1 )
860 CALL dgebr2d( ictxt, scope, top, 1, 1, alpha, 1,
864 ELSE IF( lsame( scope,
'A' ) )
THEN
865 IF( ( myrow.EQ.iarow ).AND.( mycol.EQ.iacol ) )
THEN
866 ioffa = iia+(jja-1)*desca( lld_ )
867 CALL dgebs2d( ictxt, scope, top, 1, 1, a( ioffa ), 1 )
870 CALL dgebr2d( ictxt, scope, top, 1, 1, alpha, 1,
874 IF( myrow.EQ.iarow .AND. mycol.EQ.iacol )
875 alpha = a( iia+(jja-1)*desca( lld_ ) )
884 SUBROUTINE pdelset( A, IA, JA, DESCA, ALPHA )
885 use descriptor_mod, desca_x=>desca
992 INTEGER IACOL, IAROW, IIA, JJA
995 EXTERNAL blacs_gridinfo, infog2l
1001 CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
1003 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
1006 IF( myrow.EQ.iarow .AND. mycol.EQ.iacol )
1007 a( iia+(jja-1)*desca( lld_ ) ) = alpha