From 54c98ac60817752d8f0191d86684a304973ac97c Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 29 Nov 2019 15:49:52 +0000 Subject: [PATCH] New fnd_owner Merged new fnd_owner algorithm from development/cdasb --- base/internals/Makefile | 6 +- base/internals/psb_indx_map_fnd_owner.F90 | 166 +-------- base/internals/psi_a2a_fnd_owner.F90 | 121 ++++++ base/internals/psi_adjcncy_fnd_owner.F90 | 384 ++++++++++++++++++++ base/internals/psi_fnd_owner.F90 | 2 +- base/internals/psi_graph_fnd_owner.F90 | 319 ++++++++++++++++ base/internals/psi_symm_dep_list.F90 | 162 +++++++++ base/modules/desc/psb_desc_mod.F90 | 59 ++- base/modules/desc/psb_gen_block_map_mod.f90 | 2 +- base/modules/desc/psb_glist_map_mod.f90 | 2 +- base/modules/desc/psb_indx_map_mod.f90 | 100 ++++- base/modules/desc/psb_repl_map_mod.f90 | 2 +- 12 files changed, 1153 insertions(+), 172 deletions(-) create mode 100644 base/internals/psi_a2a_fnd_owner.F90 create mode 100644 base/internals/psi_adjcncy_fnd_owner.F90 create mode 100644 base/internals/psi_graph_fnd_owner.F90 create mode 100644 base/internals/psi_symm_dep_list.F90 diff --git a/base/internals/Makefile b/base/internals/Makefile index 471fc383..f1f8d61b 100644 --- a/base/internals/Makefile +++ b/base/internals/Makefile @@ -2,11 +2,11 @@ include ../../Make.inc FOBJS = psi_compute_size.o psi_crea_bnd_elem.o psi_crea_index.o \ psi_crea_ovr_elem.o psi_bld_tmpovrl.o psi_dl_check.o \ - psi_bld_tmphalo.o psi_sort_dl.o \ + psi_bld_tmphalo.o psi_sort_dl.o psb_indx_map_fnd_owner.o \ psi_desc_impl.o psi_list_search.o psi_srtlist.o -MPFOBJS = psi_desc_index.o psi_extrct_dl.o \ - psi_fnd_owner.o psb_indx_map_fnd_owner.o +MPFOBJS = psi_desc_index.o psi_extrct_dl.o psi_fnd_owner.o psi_a2a_fnd_owner.o \ + psi_graph_fnd_owner.o psi_adjcncy_fnd_owner.o psi_symm_dep_list.o LIBDIR=.. INCDIR=.. diff --git a/base/internals/psb_indx_map_fnd_owner.F90 b/base/internals/psb_indx_map_fnd_owner.F90 index fa439062..3e2cb576 100644 --- a/base/internals/psb_indx_map_fnd_owner.F90 +++ b/base/internals/psb_indx_map_fnd_owner.F90 @@ -62,14 +62,11 @@ subroutine psb_indx_map_fnd_owner(idx,iprc,idxmap,info) #endif integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - class(psb_indx_map), intent(in) :: idxmap + class(psb_indx_map), intent(inout) :: idxmap integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), allocatable :: helem(:),hproc(:),& - & answers(:,:),idxsrch(:,:), hhidx(:) - integer(psb_mpik_), allocatable :: hsz(:),hidx(:), & - & sdsz(:),sdidx(:), rvsz(:), rvidx(:) + integer(psb_ipk_), allocatable :: hhidx(:) integer(psb_mpik_) :: icomm, minfo, iictxt integer(psb_ipk_) :: i,n_row,n_col,err_act,ih,hsize,ip,isz,k,j,& & last_ih, last_j, nv, mglob @@ -139,166 +136,11 @@ subroutine psb_indx_map_fnd_owner(idx,iprc,idxmap,info) else - ! - ! The basic idea is very simple. - ! First we collect (to all) all the requests. - Allocate(hidx(np+1),hsz(np),& - & sdsz(0:np-1),sdidx(0:np-1),& - & rvsz(0:np-1),rvidx(0:np-1),& - & stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - hsz = 0 - hsz(me+1) = nv - call psb_amx(iictxt,hsz) - hidx(1) = 0 - do i=1, np - hidx(i+1) = hidx(i) + hsz(i) - end do - hsize = hidx(np+1) - Allocate(helem(hsize),hproc(hsize),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - if (gettime) then - t3 = psb_wtime() - end if - - call mpi_allgatherv(idx,hsz(me+1),psb_mpi_ipk_integer,& - & hproc,hsz,hidx,psb_mpi_ipk_integer,& - & icomm,minfo) - if (gettime) then - tamx = psb_wtime() - t3 - end if - - ! Second, we figure out locally whether we own the indices (whoever is - ! asking for them). - if (gettime) then - t3 = psb_wtime() - end if - call idxmap%g2l(hproc(1:hsize),helem(1:hsize),info,owned=.true.) - if (gettime) then - tidx = psb_wtime()-t3 - end if - if (info == psb_err_iarray_outside_bounds_) info = psb_success_ - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_idx_cnv') - goto 9999 - end if + call psi_graph_fnd_owner(idx,iprc,idxmap,info) - ! Third: we build the answers for those indices we own, - ! with a section for each process asking. - hidx = hidx +1 - j = 0 - do ip = 0, np-1 - sdidx(ip) = j - sdsz(ip) = 0 - do i=hidx(ip+1), hidx(ip+1+1)-1 - if ((0 < helem(i)).and. (helem(i) <= n_row)) then - j = j + 1 - hproc(j) = hproc(i) - sdsz(ip) = sdsz(ip) + 1 - end if - end do - end do - - if (gettime) then - t3 = psb_wtime() - end if - - ! Collect all the answers with alltoallv (need sizes) - call mpi_alltoall(sdsz,1,psb_mpi_def_integer,& - & rvsz,1,psb_mpi_def_integer,icomm,minfo) - - isz = sum(rvsz) - - allocate(answers(isz,2),idxsrch(nv,2),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - j = 0 - do ip=0, np-1 - rvidx(ip) = j - j = j + rvsz(ip) - end do - call mpi_alltoallv(hproc,sdsz,sdidx,psb_mpi_ipk_integer,& - & answers(:,1),rvsz,rvidx,psb_mpi_ipk_integer,& - & icomm,minfo) - if (gettime) then - tamx = psb_wtime() - t3 + tamx - end if - j = 1 - do ip = 0,np-1 - do k=1,rvsz(ip) - answers(j,2) = ip - j = j + 1 - end do - end do - ! Sort the answers and the requests, so we can - ! match them efficiently - call psb_msort(answers(:,1),ix=answers(:,2),& - & flag=psb_sort_keep_idx_) - idxsrch(1:nv,1) = idx(1:nv) - call psb_msort(idxsrch(1:nv,1),ix=idxsrch(1:nv,2)) - - ! Now extract the answers for our local query - last_ih = -1 - last_j = -1 - j = 1 - do i=1, nv - ih = idxsrch(i,1) - if (ih == last_ih) then - iprc(idxsrch(i,2)) = answers(last_j,2) - else - - do - if (j > size(answers,1)) then - ! Last resort attempt. - j = psb_ibsrch(ih,size(answers,1,kind=psb_ipk_),answers(:,1)) - if (j == -1) then - write(psb_err_unit,*) me,'psi_fnd_owner: searching for ',ih, & - & 'not found : ',size(answers,1),':',answers(:,1) - info = psb_err_internal_error_ - call psb_errpush(psb_err_internal_error_,name,a_err='out bounds srch ih') - goto 9999 - end if - end if - if (answers(j,1) == ih) exit - if (answers(j,1) > ih) then - k = j - j = psb_ibsrch(ih,k,answers(1:k,1)) - if (j == -1) then - write(psb_err_unit,*) me,'psi_fnd_owner: searching for ',ih, & - & 'not found : ',size(answers,1),':',answers(:,1) - info = psb_err_internal_error_ - call psb_errpush(psb_err_internal_error_,name,a_err='out bounds srch ih') - goto 9999 - end if - end if - - j = j + 1 - end do - ! Note that the answers here are given in order - ! of sending process, so we are implicitly getting - ! the max process index in case of overlap. - last_ih = ih - do - last_j = j - iprc(idxsrch(i,2)) = answers(j,2) - j = j + 1 - if (j > size(answers,1)) exit - if (answers(j,1) /= ih) exit - end do - end if - end do end if + if (gettime) then call psb_barrier(ictxt) t1 = psb_wtime() diff --git a/base/internals/psi_a2a_fnd_owner.F90 b/base/internals/psi_a2a_fnd_owner.F90 new file mode 100644 index 00000000..e755bf82 --- /dev/null +++ b/base/internals/psi_a2a_fnd_owner.F90 @@ -0,0 +1,121 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! File: psi_fnd_owner.f90 +! +! Subroutine: psi_fnd_owner +! Figure out who owns global indices. +! +! Arguments: +! nv - integer Number of indices required on the calling +! process +! idx(:) - integer Required indices on the calling process. +! Note: the indices should be unique! +! iprc(:) - integer(psb_ipk_), allocatable Output: process identifiers for the corresponding +! indices +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. return code. +! +subroutine psi_a2a_fnd_owner(idx,iprc,idxmap,info) + use psb_serial_mod + use psb_const_mod + use psb_error_mod + use psb_penv_mod + use psb_realloc_mod + use psb_indx_map_mod, psb_protect_name => psi_a2a_fnd_owner +#ifdef MPI_MOD + use mpi +#endif + + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), allocatable, intent(out) :: iprc(:) + class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), intent(out) :: info + + + integer(psb_ipk_), allocatable :: tmpadj(:) + integer(psb_mpik_) :: icomm, minfo, iictxt + integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,isz,j, k,& + & last_ih, last_j, nv + integer(psb_ipk_) :: mglob, ih + integer(psb_ipk_) :: ictxt,np,me, nresp + logical, parameter :: gettime=.false. + real(psb_dpk_) :: t0, t1, t2, t3, t4, tamx, tidx + character(len=20) :: name + + info = psb_success_ + name = 'psi_a2a_fnd_owner' + call psb_erractionsave(err_act) + + ictxt = idxmap%get_ctxt() + icomm = idxmap%get_mpic() + mglob = idxmap%get_gr() + n_row = idxmap%get_lr() + n_col = idxmap%get_lc() + iictxt = ictxt + + call psb_info(ictxt, me, np) + + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + if (.not.(idxmap%is_valid())) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid idxmap') + goto 9999 + end if + + + ! + ! Reuse the other version by tricking it with an adjcncy list + ! that contains everybody but ME. + ! + nv = size(idx) + call psb_realloc(np-1,tmpadj,info) + tmpadj(1:me) = [(i,i=0,me-1)] + tmpadj(me+1:np-1) = [(i,i=me+1,np-1)] + call psi_adjcncy_fnd_owner(idx,iprc,tmpadj,idxmap,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psi_a2a_fnd_owner diff --git a/base/internals/psi_adjcncy_fnd_owner.F90 b/base/internals/psi_adjcncy_fnd_owner.F90 new file mode 100644 index 00000000..ad66904a --- /dev/null +++ b/base/internals/psi_adjcncy_fnd_owner.F90 @@ -0,0 +1,384 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! File: psi_fnd_owner.f90 +! +! Subroutine: psi_fnd_owner +! Figure out who owns global indices. +! +! Arguments: +! +subroutine psi_adjcncy_fnd_owner(idx,iprc,adj,idxmap,info) + use psb_serial_mod + use psb_const_mod + use psb_error_mod + use psb_penv_mod + use psb_realloc_mod + use psb_indx_map_mod, psb_protect_name => psi_adjcncy_fnd_owner +#ifdef MPI_MOD + use mpi +#endif + + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), allocatable, intent(out) :: iprc(:) + integer(psb_ipk_), intent(in) :: adj(:) + class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), intent(out) :: info + + + integer(psb_ipk_), allocatable :: rmtidx(:) + integer(psb_ipk_), allocatable :: tproc(:), lclidx(:) + integer(psb_mpik_), allocatable :: hsz(:),hidx(:), sdidx(:), rvidx(:),& + & sdsz(:), rvsz(:), sdhd(:), rvhd(:), p2pstat(:,:) + integer(psb_mpik_) :: prc, p2ptag, iret + integer(psb_mpik_) :: icomm, minfo, iictxt + integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,isz,j, k,& + & last_ih, last_j, nidx, nrecv, nadj + integer(psb_ipk_) :: mglob, ih + integer(psb_ipk_) :: ictxt,np,me + logical, parameter :: gettime=.false., new_impl=.true. + logical, parameter :: a2av_impl=.true., debug=.false. + real(psb_dpk_) :: t0, t1, t2, t3, t4, tamx, tidx + character(len=20) :: name + + info = psb_success_ + name = 'psi_adjcncy_fnd_owner' + call psb_erractionsave(err_act) + + ictxt = idxmap%get_ctxt() + icomm = idxmap%get_mpic() + mglob = idxmap%get_gr() + n_row = idxmap%get_lr() + n_col = idxmap%get_lc() + iictxt = ictxt + + call psb_info(ictxt, me, np) + + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + if (.not.(idxmap%is_valid())) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid idxmap') + goto 9999 + end if + + if (gettime) then + t0 = psb_wtime() + end if + + nadj = size(adj) + nidx = size(idx) + call psb_realloc(nidx,iprc,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_realloc') + goto 9999 + end if + iprc = -1 + ! write(0,*) me,name,' Going through ',nidx,nadj + + if (a2av_impl) then + ! + ! First simple minded version with auxiliary arrays + ! dimensioned on NP. + ! Do the exchange with an alltoallv + ! + ! + Allocate(hidx(0:np),hsz(np),sdsz(0:np-1),rvsz(0:np-1), & + & sdidx(0:np),rvidx(0:np),stat=info) + ! + ! Same send buffer for everybody + ! + sdidx(:) = 0 + ! + ! First, send sizes according to adjcncy list + ! + sdsz = 0 + do j=1, nadj + sdsz(adj(j)) = nidx + end do + !write(0,*)me,' Check on sizes into a2a:',adj(:),nadj,':',sdsz(:) + + call mpi_alltoall(sdsz,1,psb_mpi_def_integer,& + & rvsz,1,psb_mpi_def_integer,icomm,minfo) + rvidx(0) = 0 + do i=0, np-1 + rvidx(i+1) = rvidx(i) + rvsz(i) + end do + hsize = rvidx(np) + ! write(0,*)me,' Check on sizes from a2a:',hsize,rvsz(:) + ! + ! Second, allocate buffers and exchange data + ! + Allocate(rmtidx(hsize),lclidx(max(hsize,nidx*nadj)),& + & tproc(max(hsize,nidx)),stat=info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + call mpi_alltoallv(idx,sdsz,sdidx,psb_mpi_ipk_integer,& + & rmtidx,rvsz,rvidx,psb_mpi_ipk_integer,icomm,iret) + + ! + ! Third, compute local answers + ! + call idxmap%g2l(rmtidx(1:hsize),lclidx(1:hsize),info,owned=.true.) + do i=1, hsize + tproc(i) = -1 + if ((0 < lclidx(i)).and. (lclidx(i) <= n_row)) tproc(i) = me + end do + + ! + ! Fourth, exchange the answers + ! + ! Adjust sdidx for reuse in receiving lclidx array + do i=0,np-1 + sdidx(i+1) = sdidx(i) + sdsz(i) + end do + call mpi_alltoallv(tproc,rvsz,rvidx,psb_mpi_ipk_integer,& + & lclidx,sdsz,sdidx,psb_mpi_ipk_integer,icomm,iret) + + do i=0, np-1 + if (sdsz(i)>0) then + ! Must be nidx == sdsz(i) + iprc(1:nidx) = max(iprc(1:nidx), lclidx(sdidx(i)+1:sdidx(i)+sdsz(i))) + end if + end do + + if (debug) write(0,*) me,' End of adjcncy_fnd ',iprc(1:nidx) + else + if (new_impl) then + ! + ! First simple minded version with auxiliary arrays + ! dimensioned on NP. + ! Could it be improved with a loop based on the maximum length + ! of adj(:) ??? + ! + Allocate(hidx(0:np),hsz(np),sdsz(0:np-1),rvsz(0:np-1),& + & sdhd(0:np-1), rvhd(0:np-1), p2pstat(mpi_status_size,0:np-1),& + & stat=info) + sdhd(:) = mpi_request_null + rvhd(:) = mpi_request_null + ! + ! First, send sizes according to adjcncy list + ! + sdsz = 0 + do j=1, nadj + sdsz(adj(j)) = nidx + end do + !write(0,*)me,' Check on sizes into a2a:',adj(:),nadj,':',sdsz(:) + + call mpi_alltoall(sdsz,1,psb_mpi_def_integer,& + & rvsz,1,psb_mpi_def_integer,icomm,minfo) + hidx(0) = 0 + do i=0, np-1 + hidx(i+1) = hidx(i) + rvsz(i) + end do + hsize = hidx(np) + ! write(0,*)me,' Check on sizes from a2a:',hsize,rvsz(:) + ! + ! Second, allocate buffers and exchange data + ! + Allocate(rmtidx(hsize),lclidx(max(hsize,nidx*nadj)),tproc(max(hsize,nidx)),stat=info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + do i = 0, np-1 + if (rvsz(i)>0) then + ! write(0,*) me, ' First receive from ',i,rvsz(i) + call psb_get_rank(prc,ictxt,i) + p2ptag = psb_int_swap_tag + !write(0,*) me, ' Posting first receive from ',i,rvsz(i),prc + call mpi_irecv(rmtidx(hidx(i)+1),rvsz(i),& + & psb_mpi_ipk_integer,prc,& + & p2ptag, icomm,rvhd(i),iret) + end if + end do + do j=1, nadj + if (nidx > 0) then + !call psb_snd(ictxt,idx(1:nidx),adj(j)) + call psb_get_rank(prc ,ictxt,adj(j)) + p2ptag = psb_int_swap_tag + !write(0,*) me, ' First send to ',adj(j),nidx, prc + call mpi_send(idx,nidx,& + & psb_mpi_ipk_integer,prc,& + & p2ptag, icomm,iret) + end if + end do +!!$ do i = 0, np-1 +!!$ if (rvsz(i)>0) then +!!$ ! write(0,*) me, ' First receive from ',i,rvsz(i) +!!$ call psb_rcv(ictxt,rmtidx(hidx(i)+1:hidx(i)+rvsz(i)),i) +!!$ end if +!!$ end do + call mpi_waitall(np,rvhd,p2pstat,iret) + + ! + ! Third, compute local answers + ! + call idxmap%g2l(rmtidx(1:hsize),lclidx(1:hsize),info,owned=.true.) + do i=1, hsize + tproc(i) = -1 + if ((0 < lclidx(i)).and. (lclidx(i) <= n_row)) tproc(i) = me + end do + ! + ! At this point we can reuse lclidx to receive messages + ! + rvhd(:) = mpi_request_null + do j=1, nadj + !write(0,*) me, ' First send to ',adj(j),nidx + if (nidx > 0) then + !call psb_snd(ictxt,idx(1:nidx),adj(j)) + call psb_get_rank(prc,ictxt,adj(j)) + p2ptag = psb_int_swap_tag + !write(0,*) me, ' Posting second receive from ',adj(j),nidx, prc + call mpi_irecv(lclidx((j-1)*nidx+1),nidx, & + & psb_mpi_ipk_integer,prc,& + & p2ptag, icomm,rvhd(j),iret) + end if + end do + + ! + ! Fourth, send data back; + ! + do i = 0, np-1 + if (rvsz(i)>0) then + !call psb_snd(ictxt,tproc(hidx(i)+1:hidx(i)+rvsz(i)),i) + call psb_get_rank(prc,ictxt,i) + p2ptag = psb_int_swap_tag + !write(0,*) me, ' Second send to ',i,rvsz(i), prc + call mpi_send(tproc(hidx(i)+1),rvsz(i),& + & psb_mpi_ipk_integer,prc,& + & p2ptag, icomm,iret) + end if + end do + ! + ! Fifth: receive and combine. MAX works because default + ! answer is -1. + ! + call mpi_waitall(np,rvhd,p2pstat,iret) + do j = 1, nadj + !write(0,*) me, ' Second receive from ',adj(j), nidx + !if (nidx > 0) call psb_rcv(ictxt,tproc(1:nidx),adj(j)) + iprc(1:nidx) = max(iprc(1:nidx), lclidx((j-1)*nidx+1:(j-1)*nidx+nidx)) + end do + if (debug) write(0,*) me,' End of adjcncy_fnd ',iprc(1:nidx) + + else + + Allocate(hidx(0:np),hsz(np),& + & sdsz(0:np-1),rvsz(0:np-1),stat=info) + ! + ! First, send sizes according to adjcncy list + ! + sdsz = 0 + do j=1, nadj + sdsz(adj(j)) = nidx + end do + !write(0,*)me,' Check on sizes into a2a:',adj(:),nadj,':',sdsz(:) + + call mpi_alltoall(sdsz,1,psb_mpi_def_integer,& + & rvsz,1,psb_mpi_def_integer,icomm,minfo) + hidx(0) = 0 + do i=0, np-1 + hidx(i+1) = hidx(i) + rvsz(i) + end do + hsize = hidx(np) + ! write(0,*)me,' Check on sizes from a2a:',hsize,rvsz(:) + ! + ! Second, allocate buffers and exchange data + ! + Allocate(rmtidx(hsize),lclidx(hsize),tproc(max(hsize,nidx)),stat=info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + do j=1, nadj + !write(0,*) me, ' First send to ',adj(j),nidx + if (nidx > 0) call psb_snd(ictxt,idx(1:nidx),adj(j)) + end do + do i = 0, np-1 + if (rvsz(i)>0) then + ! write(0,*) me, ' First receive from ',i,rvsz(i) + call psb_rcv(ictxt,rmtidx(hidx(i)+1:hidx(i)+rvsz(i)),i) + end if + end do + + ! + ! Third, compute local answers + ! + call idxmap%g2l(rmtidx(1:hsize),lclidx(1:hsize),info,owned=.true.) + do i=1, hsize + tproc(i) = -1 + if ((0 < lclidx(i)).and. (lclidx(i) <= n_row)) tproc(i) = me + end do + + ! + ! Fourth, send data back; + ! + do i = 0, np-1 + if (rvsz(i)>0) then + !write(0,*) me, ' Second send to ',i,rvsz(i) + call psb_snd(ictxt,tproc(hidx(i)+1:hidx(i)+rvsz(i)),i) + end if + end do + ! + ! Fifth: receive and combine. MAX works because default + ! answer is -1. Reuse tproc + ! + do j = 1, nadj + !write(0,*) me, ' Second receive from ',adj(j), nidx + if (nidx > 0) call psb_rcv(ictxt,tproc(1:nidx),adj(j)) + iprc(1:nidx) = max(iprc(1:nidx), tproc(1:nidx)) + end do + end if + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psi_adjcncy_fnd_owner diff --git a/base/internals/psi_fnd_owner.F90 b/base/internals/psi_fnd_owner.F90 index 3fe9d3b3..ea0bf46d 100644 --- a/base/internals/psi_fnd_owner.F90 +++ b/base/internals/psi_fnd_owner.F90 @@ -65,7 +65,7 @@ subroutine psi_fnd_owner(nv,idx,iprc,desc,info) integer(psb_ipk_), intent(in) :: nv integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - type(psb_desc_type), intent(in) :: desc + type(psb_desc_type), intent(inout) :: desc integer(psb_ipk_), intent(out) :: info diff --git a/base/internals/psi_graph_fnd_owner.F90 b/base/internals/psi_graph_fnd_owner.F90 new file mode 100644 index 00000000..3b5e85bf --- /dev/null +++ b/base/internals/psi_graph_fnd_owner.F90 @@ -0,0 +1,319 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! File: psi_fnd_owner.f90 +! +! Subroutine: psi_fnd_owner +! Figure out who owns global indices. +! +! Arguments: +! nv - integer Number of indices required on the calling +! process +! idx(:) - integer Required indices on the calling process. +! Note: the indices should be unique! +! iprc(:) - integer(psb_ipk_), allocatable Output: process identifiers for the corresponding +! indices +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. return code. +! +subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info) + use psb_serial_mod + use psb_const_mod + use psb_error_mod + use psb_penv_mod + use psb_realloc_mod + use psb_desc_mod, psb_protect_name => psi_graph_fnd_owner +#ifdef MPI_MOD + use mpi +#endif + + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), allocatable, intent(out) :: iprc(:) + class(psb_indx_map), intent(inout) :: idxmap + integer(psb_ipk_), intent(out) :: info + + + integer(psb_ipk_), allocatable :: tidx(:) + integer(psb_ipk_), allocatable :: tprc(:), tsmpl(:), ladj(:) + integer(psb_mpik_) :: icomm, minfo, iictxt + integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,isz,j,ipnt, nsampl_out,& + & last_ih, last_j, nv, n_answers, n_rest, nsampl_in, locr_max, & + & nrest_max, nadj, maxspace, mxnsin + integer(psb_ipk_) :: mglob, ih + integer(psb_ipk_) :: ictxt,np,me, nresp + integer(psb_ipk_), parameter :: nt=4 + integer(psb_ipk_) :: tmpv(2) + logical, parameter :: gettime=.false., trace=.false. + real(psb_dpk_) :: t0, t1, t2, t3, t4 + character(len=20) :: name + + info = psb_success_ + name = 'psi_graph_fnd_owner' + call psb_erractionsave(err_act) + + ictxt = idxmap%get_ctxt() + icomm = idxmap%get_mpic() + mglob = idxmap%get_gr() + n_row = idxmap%get_lr() + n_col = idxmap%get_lc() + iictxt = ictxt + + call psb_info(ictxt, me, np) + + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + if (.not.(idxmap%is_valid())) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid idxmap') + goto 9999 + end if + + ! + ! Choice of maxspace should be adjusted to account for a default + ! "sensible" size and/or a user-specified value + ! + tmpv(1) = n_row + tmpv(2) = psb_cd_get_maxspace() + call psb_max(ictxt,tmpv) + locr_max = tmpv(1) + maxspace = min(nt*locr_max,tmpv(2)) + maxspace = max(maxspace,np) + if (trace.and.(me == 0)) write(0,*) ' Through graph_fnd_owner with maxspace:',maxspace + ! + ! + ! + nv = size(idx) + call psb_realloc(nv,iprc,info) + if (info == psb_success_) call psb_realloc(nv,tidx,info) + if (info == psb_success_) call psb_realloc(nv,tprc,info) + if (info == psb_success_) call psb_realloc(nv,tsmpl,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_realloc') + goto 9999 + end if + iprc(:) = -1 + n_answers = 0 + ! + ! Start from the adjacncy list + ! + call psb_safe_ab_cpy(idxmap%p_adjcncy,ladj,info) + nadj = psb_size(ladj) + ! This makes ladj allocated with size 0 just in case + call psb_realloc(nadj,ladj,info) + n_rest = nv - n_answers + nrest_max = n_rest + + tmpv(1) = nadj + tmpv(2) = nrest_max + call psb_max(ictxt,tmpv) + if ((tmpv(1) > 0).and.(tmpv(2) >0)) then + ! + ! Do a preliminary run on the user-defined adjacency lists + ! + if (trace.and.(me == 0)) write(0,*) ' Initial sweep on user-defined topology' + nsampl_in = min(n_rest,max(1,(maxspace+max(1,nadj)-1))/(max(1,nadj))) + call psi_adj_fnd_sweep(idx,iprc,ladj,idxmap,nsampl_in,n_answers) + call idxmap%xtnd_p_adjcncy(ladj) + n_rest = nv - n_answers + nrest_max = n_rest + call psb_max(ictxt,nrest_max) + if (trace.and.(me == 0)) write(0,*) ' After initial sweep:',nrest_max + end if + + + fnd_owner_loop: do while (nrest_max>0) + ! + ! The basic idea of this loop is to alternate between + ! searching through all processes and searching + ! in the neighbourood. + ! + ! 1. Select a sample such that the total size is <= maxspace + ! sample query is then sent to all processes + ! + ! if (trace.and.(me == 0)) write(0,*) 'Looping in graph_fnd_owner: ', nrest_max + nsampl_in = min(n_rest,max(1,(maxspace+np-1)/np)) + ! + ! Choose a sample, should it be done in this simplistic way? + ! Note: nsampl_in is a hint, not an absolute, hence nsampl_out + ! + ipnt = 1 +!!$ write(0,*) me,' Into first sampling ',nsampl_in + call psi_get_sample(ipnt, idx,iprc,tidx,tsmpl,nsampl_in,nsampl_out) + nsampl_in = min(nsampl_out,nsampl_in) +!!$ write(0,*) me,' From first sampling ',nsampl_in + ! + ! 2. Do a search on all processes; this is supposed to find + ! the owning process for all inputs; + ! + call psi_a2a_fnd_owner(tidx(1:nsampl_in),tprc,idxmap,info) + call psi_cpy_out(iprc,tprc,tsmpl,nsampl_in,nsampl_out) + if (nsampl_out /= nsampl_in) then + write(0,*) me,'Warning: indices not found by a2a_fnd_owner ',nsampl_out,nsampl_in + end if + n_answers = n_answers + nsampl_out + n_rest = nv - n_answers + ! + ! 3. Extract the resulting adjacency list and add it to the + ! indxmap; + ! + ladj = tprc(1:nsampl_in) + call psb_msort_unique(ladj,nadj) + call psb_realloc(nadj,ladj,info) + + ! + ! 4. Extract again a sample and do a neighbourhood search + ! so that the total size is <= maxspace + ! (will not be exact since nadj varies with process) + ! Need to set up a proper loop here to have a complete + ! sweep over the input vector. Done inside adj_fnd_sweep. + ! +!!$ write(0,*) me,' After a2a ',n_rest + nsampl_in = min(n_rest,max(1,(maxspace+max(1,nadj)-1))/(max(1,nadj))) + mxnsin = nsampl_in + call psb_max(ictxt,mxnsin) +!!$ write(0,*) me, ' mxnsin ',mxnsin + if (mxnsin>0) call psi_adj_fnd_sweep(idx,iprc,ladj,idxmap,nsampl_in,n_answers) + call idxmap%xtnd_p_adjcncy(ladj) + + n_rest = nv - n_answers + nrest_max = n_rest + call psb_max(ictxt,nrest_max) + if (trace.and.(me == 0)) write(0,*) ' fnd_owner_loop remaining:',nrest_max + end do fnd_owner_loop + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +contains + + subroutine psi_get_sample(ipntidx,idx,iprc,tidx,tsmpl,ns_in,ns_out) + implicit none + integer(psb_ipk_), intent(inout) :: ipntidx + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), intent(in) :: ns_in, iprc(:) + integer(psb_ipk_), intent(out) :: tidx(:) + integer(psb_ipk_), intent(out) :: tsmpl(:), ns_out + ! + integer(psb_ipk_) :: nv, ns + + nv = size(idx) + ! + ! Choose a sample, should it be done in this simplistic way? + ! + ns = ns_in + ! + ! ns_in == 0 means that on the outside we figure there's + ! nothing left, but we are here because we have to synchronize. + ! Make sure we sweep through the entire vector immediately + ! + if (ns == 0) ns = nv + ns_out = 0 + + do while (ipntidx<= nv) + if (iprc(ipntidx) == -1) then + ns_out = ns_out + 1 + tsmpl(ns_out) = ipntidx + tidx(ns_out) = idx(ipntidx) + end if + ipntidx = ipntidx + 1 + if (ns_out >= ns) exit + end do + + end subroutine psi_get_sample + + subroutine psi_cpy_out(iprc,tprc,tsmpl,ns_in,ns_out) + implicit none + integer(psb_ipk_), intent(out) :: iprc(:) + integer(psb_ipk_), intent(in) :: ns_in + integer(psb_ipk_), intent(in) :: tprc(:), tsmpl(:) + integer(psb_ipk_), intent(out) :: ns_out + + integer(psb_ipk_) :: j + + ns_out = 0 + do j=1, ns_in + if (tprc(j) /= -1) then + ns_out = ns_out + 1 + iprc(tsmpl(j)) = tprc(j) + end if + end do + end subroutine psi_cpy_out + + subroutine psi_adj_fnd_sweep(idx,iprc,adj,idxmap,n_samples,n_answers) + implicit none + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), intent(in) :: n_samples + integer(psb_ipk_), intent(inout) :: iprc(:), n_answers + integer(psb_ipk_), intent(in) :: adj(:) + class(psb_indx_map), intent(inout) :: idxmap + ! + integer(psb_ipk_) :: ipnt, ns_in, ns_out, n_rem, ictxt, me, np, isw + integer(psb_ipk_), allocatable :: tidx(:) + integer(psb_ipk_), allocatable :: tsmpl(:) + + ictxt = idxmap%get_ctxt() + call psb_info(ictxt,me,np) + call psb_realloc(n_samples,tidx,info) + call psb_realloc(n_samples,tsmpl,info) + ipnt = 1 + isw = 1 + do + !write(0,*) me,' Into sampling ',n_samples + call psi_get_sample(ipnt, idx,iprc,tidx,tsmpl,n_samples,ns_out) + ns_in = min(n_samples,ns_out) + !write(0,*) me,' From second sampling ',ns_out + call psi_adjcncy_fnd_owner(tidx(1:ns_in),tprc,ladj,idxmap,info) + call psi_cpy_out(iprc,tprc,tsmpl,ns_in,ns_out) + !write(0,*) me,' Sweep ',isw,' answers:',ns_out + n_answers = n_answers + ns_out + n_rem = size(idx)-ipnt + call psb_max(ictxt,n_rem) + !write(0,*) me,' Sweep ',isw,n_rem, ipnt, n_samples + if (n_rem <= 0) exit + isw = isw + 1 + end do + + + end subroutine psi_adj_fnd_sweep + +end subroutine psi_graph_fnd_owner diff --git a/base/internals/psi_symm_dep_list.F90 b/base/internals/psi_symm_dep_list.F90 new file mode 100644 index 00000000..4cb1001d --- /dev/null +++ b/base/internals/psi_symm_dep_list.F90 @@ -0,0 +1,162 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! File: psi_fnd_owner.f90 +! +! Subroutine: psi_fnd_owner +! Figure out who owns global indices. +! +! Arguments: +! +subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) + use psb_serial_mod + use psb_const_mod + use psb_error_mod + use psb_penv_mod + use psb_realloc_mod + use psb_indx_map_mod, psb_protect_name => psi_symm_dep_list +#ifdef MPI_MOD + use mpi +#endif + + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + integer(psb_mpik_), intent(inout) :: rvsz(0:) + integer(psb_ipk_), allocatable, intent(inout) :: adj(:) + class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: flag + + ! + integer(psb_ipk_), allocatable :: ladj(:) + integer(psb_mpik_), allocatable :: sdsz(:) + integer(psb_mpik_) :: icomm, minfo, iictxt + integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,& + & last_ih, last_j, nidx, nrecv, nadj, flag_ + integer(psb_ipk_) :: mglob, ih + integer(psb_ipk_) :: ictxt,np,me + character(len=20) :: name + + info = psb_success_ + name = 'psi_symm_dep_list' + call psb_erractionsave(err_act) + + ictxt = idxmap%get_ctxt() + icomm = idxmap%get_mpic() + mglob = idxmap%get_gr() + n_row = idxmap%get_lr() + n_col = idxmap%get_lc() + iictxt = ictxt + + call psb_info(ictxt, me, np) + + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + if (.not.(idxmap%is_valid())) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid idxmap') + goto 9999 + end if + + + nadj = size(adj) + if (size(rvsz) 0) nrecv = nrecv + 1 + end do + + ! + ! Now fix adj to be symmetric + ! + call psb_realloc(nadj+nrecv,ladj,info) + ladj(1:nadj) = adj(1:nadj) + do ip=0, np-1 + if (rvsz(ip)>0) then + nadj = nadj + 1 + ladj(nadj+1:nadj+nrecv) = ip + end if + end do + call psb_msort_unique(ladj,nadj) + call psb_realloc(nadj,adj,info) + adj(1:nadj) = ladj(1:nadj) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psi_symm_dep_list diff --git a/base/modules/desc/psb_desc_mod.F90 b/base/modules/desc/psb_desc_mod.F90 index cfec78b3..ea395615 100644 --- a/base/modules/desc/psb_desc_mod.F90 +++ b/base/modules/desc/psb_desc_mod.F90 @@ -233,6 +233,9 @@ module psb_desc_mod procedure, pass(desc) :: get_global_rows => psb_cd_get_global_rows procedure, pass(desc) :: get_global_cols => psb_cd_get_global_cols procedure, pass(desc) :: get_global_indices => psb_cd_get_global_indices + procedure, pass(desc) :: get_p_adjcncy => cd_get_p_adjcncy + procedure, pass(desc) :: set_p_adjcncy => cd_set_p_adjcncy + procedure, pass(desc) :: xtnd_p_adjcncy => cd_xtnd_p_adjcncy procedure, pass(desc) :: a_get_list => psb_cd_get_list procedure, pass(desc) :: v_get_list => psb_cd_v_get_list generic, public :: get_list => a_get_list, v_get_list @@ -278,6 +281,15 @@ module psb_desc_mod module procedure psb_cdtransfer end interface psb_move_alloc + + interface psb_cd_set_maxspace + module procedure psb_cd_set_maxspace + end interface psb_cd_set_maxspace + + interface psb_cd_get_maxspace + module procedure psb_cd_get_maxspace + end interface psb_cd_get_maxspace + private :: nullify_desc, cd_get_fmt,& & cd_l2gs1, cd_l2gs2, cd_l2gv1, cd_l2gv2, cd_g2ls1,& & cd_g2ls2, cd_g2lv1, cd_g2lv2, cd_g2ls1_ins,& @@ -285,6 +297,7 @@ module psb_desc_mod integer(psb_ipk_), private, save :: cd_large_threshold=psb_default_large_threshold + integer(psb_ipk_), private, save :: cd_maxspace = 1000*1000 contains @@ -329,6 +342,21 @@ contains val = cd_large_threshold end function psb_cd_get_large_threshold + + subroutine psb_cd_set_maxspace(ith) + implicit none + integer(psb_ipk_), intent(in) :: ith + if (ith > 0) then + cd_maxspace = ith + end if + end subroutine psb_cd_set_maxspace + + function psb_cd_get_maxspace() result(val) + implicit none + integer(psb_ipk_) :: val + val = cd_maxspace + end function psb_cd_get_maxspace + function psb_cd_is_large_size(m) result(val) use psb_penv_mod @@ -598,6 +626,35 @@ contains end function psb_cd_get_mpic + function cd_get_p_adjcncy(desc) result(val) + implicit none + integer(psb_ipk_), allocatable :: val(:) + class(psb_desc_type), intent(in) :: desc + + if (allocated(desc%indxmap)) then + val = desc%indxmap%get_p_adjcncy() + endif + + end function cd_get_p_adjcncy + + subroutine cd_set_p_adjcncy(desc,val) + implicit none + class(psb_desc_type), intent(inout) :: desc + integer(psb_ipk_), intent(in) :: val(:) + if (allocated(desc%indxmap)) then + call desc%indxmap%xtnd_p_adjcncy(val) + endif + end subroutine cd_set_p_adjcncy + + subroutine cd_xtnd_p_adjcncy(desc,val) + implicit none + class(psb_desc_type), intent(inout) :: desc + integer(psb_ipk_), intent(in) :: val(:) + if (allocated(desc%indxmap)) then + call desc%indxmap%xtnd_p_adjcncy(val) + endif + end subroutine cd_xtnd_p_adjcncy + subroutine psb_cd_set_ovl_asb(desc,info) ! ! Change state of a descriptor into ovl_build. @@ -1603,7 +1660,7 @@ contains implicit none integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - class(psb_desc_type), intent(in) :: desc + class(psb_desc_type), intent(inout) :: desc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act character(len=20) :: name='cd_fnd_owner' diff --git a/base/modules/desc/psb_gen_block_map_mod.f90 b/base/modules/desc/psb_gen_block_map_mod.f90 index 0a40d4b1..2c1ddfa4 100644 --- a/base/modules/desc/psb_gen_block_map_mod.f90 +++ b/base/modules/desc/psb_gen_block_map_mod.f90 @@ -1011,7 +1011,7 @@ contains implicit none integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - class(psb_gen_block_map), intent(in) :: idxmap + class(psb_gen_block_map), intent(inout) :: idxmap integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: ictxt, iam, np, nv, ip, i diff --git a/base/modules/desc/psb_glist_map_mod.f90 b/base/modules/desc/psb_glist_map_mod.f90 index 1e96fb17..4913702d 100644 --- a/base/modules/desc/psb_glist_map_mod.f90 +++ b/base/modules/desc/psb_glist_map_mod.f90 @@ -155,7 +155,7 @@ contains implicit none integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - class(psb_glist_map), intent(in) :: idxmap + class(psb_glist_map), intent(inout) :: idxmap integer(psb_ipk_), intent(out) :: info integer(psb_mpik_) :: ictxt, iam, np integer(psb_ipk_) :: nv, i, ngp diff --git a/base/modules/desc/psb_indx_map_mod.f90 b/base/modules/desc/psb_indx_map_mod.f90 index 96d39482..6af3b234 100644 --- a/base/modules/desc/psb_indx_map_mod.f90 +++ b/base/modules/desc/psb_indx_map_mod.f90 @@ -127,7 +127,9 @@ module psb_indx_map_mod integer(psb_ipk_), allocatable :: oracle(:,:) !> Halo owners integer(psb_ipk_), allocatable :: halo_owner(:) - + !> Adjacency list for processes + integer(psb_ipk_), allocatable :: p_adjcncy(:) + contains procedure, pass(idxmap) :: get_state => base_get_state @@ -149,6 +151,11 @@ module psb_indx_map_mod procedure, pass(idxmap) :: set_null => base_set_null procedure, nopass :: row_extendable => base_row_extendable + procedure, pass(idxmap) :: get_p_adjcncy => base_get_p_adjcncy + procedure, pass(idxmap) :: set_p_adjcncy => base_set_p_adjcncy + procedure, pass(idxmap) :: xtnd_p_adjcncy => base_xtnd_p_adjcncy + + procedure, pass(idxmap) :: set_gr => base_set_gr procedure, pass(idxmap) :: set_gc => base_set_gc procedure, pass(idxmap) :: set_lr => base_set_lr @@ -234,11 +241,60 @@ module psb_indx_map_mod implicit none integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - class(psb_indx_map), intent(in) :: idxmap + class(psb_indx_map), intent(inout) :: idxmap integer(psb_ipk_), intent(out) :: info end subroutine psb_indx_map_fnd_owner end interface + + interface + subroutine psi_a2a_fnd_owner(idx,iprc,idxmap,info) + import :: psb_indx_map, psb_ipk_ + implicit none + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), allocatable, intent(out) :: iprc(:) + class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), intent(out) :: info + end subroutine psi_a2a_fnd_owner + end interface + + interface + subroutine psi_adjcncy_fnd_owner(idx,iprc,adj,idxmap,info) + import :: psb_indx_map, psb_ipk_ + implicit none + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), allocatable, intent(out) :: iprc(:) + integer(psb_ipk_), allocatable, intent(inout) :: adj(:) + class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), intent(out) :: info + end subroutine psi_adjcncy_fnd_owner + end interface + + interface + subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info) + import :: psb_indx_map, psb_ipk_ + implicit none + integer(psb_ipk_), intent(in) :: idx(:) + integer(psb_ipk_), allocatable, intent(out) :: iprc(:) + class(psb_indx_map), intent(inout) :: idxmap + integer(psb_ipk_), intent(out) :: info + end subroutine psi_graph_fnd_owner + end interface + + integer, parameter :: psi_symm_flag_norv_ = 0 + integer, parameter :: psi_symm_flag_inrv_ = 1 + interface + subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) + import :: psb_indx_map, psb_ipk_, psb_mpik_ + implicit none + integer(psb_mpik_), intent(inout) :: rvsz(:) + integer(psb_ipk_), intent(in) :: adj(:) + class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: flag + end subroutine psi_symm_dep_list + end interface + contains @@ -301,6 +357,16 @@ contains end function base_get_lc + function base_get_p_adjcncy(idxmap) result(val) + use psb_realloc_mod + implicit none + class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), allocatable :: val(:) + integer(psb_ipk_) :: info + + call psb_safe_ab_cpy(idxmap%p_adjcncy,val,info) + + end function base_get_p_adjcncy function base_get_ctxt(idxmap) result(val) implicit none @@ -370,6 +436,36 @@ contains idxmap%local_cols = val end subroutine base_set_lc + + subroutine base_set_p_adjcncy(idxmap,val) + use psb_realloc_mod + use psb_sort_mod + implicit none + class(psb_indx_map), intent(inout) :: idxmap + integer(psb_ipk_), intent(in) :: val(:) + + call idxmap%xtnd_p_adjcncy(val) + + end subroutine base_set_p_adjcncy + + subroutine base_xtnd_p_adjcncy(idxmap,val) + use psb_realloc_mod + use psb_sort_mod + implicit none + class(psb_indx_map), intent(inout) :: idxmap + integer(psb_ipk_), intent(in) :: val(:) + integer(psb_ipk_) :: info, nv, nx + + nv = size(val) + nx = psb_size(idxmap%p_adjcncy) + call psb_realloc(nv+nx,idxmap%p_adjcncy,info) + idxmap%p_adjcncy(nx+1:nx+nv) = val(1:nv) + nx = size(idxmap%p_adjcncy) + call psb_msort_unique(idxmap%p_adjcncy,nx,dir=psb_sort_up_) + call psb_realloc(nx,idxmap%p_adjcncy,info) + + end subroutine base_xtnd_p_adjcncy + subroutine base_set_mpic(idxmap,val) implicit none class(psb_indx_map), intent(inout) :: idxmap diff --git a/base/modules/desc/psb_repl_map_mod.f90 b/base/modules/desc/psb_repl_map_mod.f90 index 016d2acc..6e48ae45 100644 --- a/base/modules/desc/psb_repl_map_mod.f90 +++ b/base/modules/desc/psb_repl_map_mod.f90 @@ -671,7 +671,7 @@ contains implicit none integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - class(psb_repl_map), intent(in) :: idxmap + class(psb_repl_map), intent(inout) :: idxmap integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: nv integer(psb_mpik_) :: ictxt, iam, np