From 2090a011db57c75d9715b37ba288d7fed1e371fa Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 6 Nov 2020 16:08:54 +0100 Subject: [PATCH] Include X_remap. --- base/modules/tools/psb_c_tools_mod.F90 | 14 ++ base/modules/tools/psb_d_tools_mod.F90 | 2 + base/modules/tools/psb_s_tools_mod.F90 | 14 ++ base/modules/tools/psb_z_tools_mod.F90 | 14 ++ base/tools/Makefile | 6 +- base/tools/psb_c_remap.F90 | 244 +++++++++++++++++++++++++ base/tools/psb_d_remap.F90 | 46 +++-- base/tools/psb_dasb.f90 | 27 +-- base/tools/psb_s_remap.F90 | 244 +++++++++++++++++++++++++ base/tools/psb_z_remap.F90 | 244 +++++++++++++++++++++++++ 10 files changed, 823 insertions(+), 32 deletions(-) create mode 100644 base/tools/psb_c_remap.F90 create mode 100644 base/tools/psb_s_remap.F90 create mode 100644 base/tools/psb_z_remap.F90 diff --git a/base/modules/tools/psb_c_tools_mod.F90 b/base/modules/tools/psb_c_tools_mod.F90 index 81e78d3a..6a7f91af 100644 --- a/base/modules/tools/psb_c_tools_mod.F90 +++ b/base/modules/tools/psb_c_tools_mod.F90 @@ -431,5 +431,19 @@ Module psb_c_tools_mod end function end interface + interface psb_remap + subroutine psb_c_remap(np_remap, desc_in, a_in, desc_out, a_out, info) + import + implicit none + !....parameters... + integer(psb_ipk_), intent(in) :: np_remap + type(psb_desc_type), intent(inout) :: desc_in + type(psb_cspmat_type), intent(inout) :: a_in + type(psb_cspmat_type), intent(out) :: a_out + type(psb_desc_type), intent(out) :: desc_out + integer(psb_ipk_), intent(out) :: info + end subroutine psb_c_remap + end interface psb_remap + end module psb_c_tools_mod diff --git a/base/modules/tools/psb_d_tools_mod.F90 b/base/modules/tools/psb_d_tools_mod.F90 index 4eca4899..da24e552 100644 --- a/base/modules/tools/psb_d_tools_mod.F90 +++ b/base/modules/tools/psb_d_tools_mod.F90 @@ -444,4 +444,6 @@ Module psb_d_tools_mod integer(psb_ipk_), intent(out) :: info end subroutine psb_d_remap end interface psb_remap + + end module psb_d_tools_mod diff --git a/base/modules/tools/psb_s_tools_mod.F90 b/base/modules/tools/psb_s_tools_mod.F90 index 2b6058da..754b6062 100644 --- a/base/modules/tools/psb_s_tools_mod.F90 +++ b/base/modules/tools/psb_s_tools_mod.F90 @@ -431,5 +431,19 @@ Module psb_s_tools_mod end function end interface + interface psb_remap + subroutine psb_s_remap(np_remap, desc_in, a_in, desc_out, a_out, info) + import + implicit none + !....parameters... + integer(psb_ipk_), intent(in) :: np_remap + type(psb_desc_type), intent(inout) :: desc_in + type(psb_sspmat_type), intent(inout) :: a_in + type(psb_sspmat_type), intent(out) :: a_out + type(psb_desc_type), intent(out) :: desc_out + integer(psb_ipk_), intent(out) :: info + end subroutine psb_s_remap + end interface psb_remap + end module psb_s_tools_mod diff --git a/base/modules/tools/psb_z_tools_mod.F90 b/base/modules/tools/psb_z_tools_mod.F90 index 09997e94..ca73ab72 100644 --- a/base/modules/tools/psb_z_tools_mod.F90 +++ b/base/modules/tools/psb_z_tools_mod.F90 @@ -431,5 +431,19 @@ Module psb_z_tools_mod end function end interface + interface psb_remap + subroutine psb_z_remap(np_remap, desc_in, a_in, desc_out, a_out, info) + import + implicit none + !....parameters... + integer(psb_ipk_), intent(in) :: np_remap + type(psb_desc_type), intent(inout) :: desc_in + type(psb_zspmat_type), intent(inout) :: a_in + type(psb_zspmat_type), intent(out) :: a_out + type(psb_desc_type), intent(out) :: desc_out + integer(psb_ipk_), intent(out) :: info + end subroutine psb_z_remap + end interface psb_remap + end module psb_z_tools_mod diff --git a/base/tools/Makefile b/base/tools/Makefile index 4a93d9b0..579f330b 100644 --- a/base/tools/Makefile +++ b/base/tools/Makefile @@ -5,7 +5,7 @@ FOBJS = psb_cdall.o psb_cdals.o psb_cdalv.o psb_cd_inloc.o psb_cdins.o psb_cdprt psb_cdcpy.o psb_cd_reinit.o psb_cd_switch_ovl_indxmap.o psb_cd_renum_block.o \ psb_dspalloc.o psb_dspasb.o psb_d_remap.o \ psb_dspfree.o psb_dspins.o psb_dsprn.o \ - psb_sspalloc.o psb_sspasb.o \ + psb_sspalloc.o psb_sspasb.o psb_s_remap.o \ psb_sspfree.o psb_sspins.o psb_ssprn.o\ psb_glob_to_loc.o psb_loc_to_glob.o\ psb_iallc.o psb_iasb.o psb_ifree.o psb_iins.o \ @@ -20,9 +20,9 @@ FOBJS = psb_cdall.o psb_cdals.o psb_cdalv.o psb_cd_inloc.o psb_cdins.o psb_cdprt psb_dallc_a.o psb_dasb_a.o psb_dfree_a.o psb_dins_a.o \ psb_callc_a.o psb_casb_a.o psb_cfree_a.o psb_cins_a.o \ psb_zallc_a.o psb_zasb_a.o psb_zfree_a.o psb_zins_a.o \ - psb_zspalloc.o psb_zspasb.o psb_zspfree.o\ + psb_zspalloc.o psb_zspasb.o psb_z_remap.o psb_zspfree.o\ psb_zspins.o psb_zsprn.o \ - psb_cspalloc.o psb_cspasb.o psb_cspfree.o\ + psb_cspalloc.o psb_cspasb.o psb_c_remap.o psb_cspfree.o\ psb_cspins.o psb_csprn.o psb_cd_set_bld.o \ psb_s_map.o psb_d_map.o psb_c_map.o psb_z_map.o \ psb_s_par_csr_spspmm.o psb_d_par_csr_spspmm.o psb_c_par_csr_spspmm.o psb_z_par_csr_spspmm.o \ diff --git a/base/tools/psb_c_remap.F90 b/base/tools/psb_c_remap.F90 new file mode 100644 index 00000000..54cc90df --- /dev/null +++ b/base/tools/psb_c_remap.F90 @@ -0,0 +1,244 @@ +! +! 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. +! +! +! +! Subroutine: psb_c_remap +! +! Arguments: +! desc_in - type(psb_desc_type). The communication descriptor to be cloned. +! desc_out - type(psb_desc_type). The output communication descriptor. +! info - integer. Return code. +subroutine psb_c_remap(np_remap, desc_in, a_in, desc_out, a_out, info) + + use psb_base_mod, psb_protect_name => psb_c_remap + + implicit none + !....parameters... + integer(psb_ipk_), intent(in) :: np_remap + type(psb_desc_type), intent(inout) :: desc_in + type(psb_cspmat_type), intent(inout) :: a_in + type(psb_cspmat_type), intent(out) :: a_out + type(psb_desc_type), intent(out) :: desc_out + integer(psb_ipk_), intent(out) :: info + + + ! locals + integer(psb_ipk_) :: np, me, ictxt, err_act + integer(psb_ipk_) :: rnp, rme, newctxt + integer(psb_ipk_) :: ipdest, ipd, id1, id2, imd, i, nsrc + integer(psb_ipk_), allocatable :: newnl(:), isrc(:), nzsrc(:), ids(:) + type(psb_lc_coo_sparse_mat) :: acoo_snd, acoo_rcv + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name + + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if (psb_get_errstatus() /= 0) return + info = psb_success_ + call psb_erractionsave(err_act) + name = 'psb_cdcpy' + + ictxt = desc_in%get_context() + + ! check on blacs grid + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),': Entered' + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + write(0,*) ' Remapping from ',np,' onto ', np_remap + + if (desc_in%get_fmt() == 'BLOCK') then + ! + ! Should we spread the processes in the new context, + ! or should we keep them close? + ! + if (.true.) then + allocate(ids(0:np_remap-1)) + if (np_remap <= np/2) then + ids(0) = 0 + do ipdest=1,np_remap -1 + ids(ipdest) = ids(ipdest-1) + np/np_remap + end do + write(0,*) ' IDS ',ids(:) + else + do ipdest = 0, np_remap-1 + ids(ipdest) = ipdest + end do + end if + call psb_init(newctxt,np=np_remap,basectxt=ictxt,ids=ids) + else + call psb_init(newctxt,np=np_remap,basectxt=ictxt) + end if + + call psb_info(newctxt,rme,rnp) + write(0,*) 'Old context: ',me,np,' New context: ',rme,rnp + call psb_bcast(ictxt,rnp) + allocate(newnl(rnp),stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + if (rnp >= np) then + write(0,*) ' No remapping on larger proc count now' + info = psb_err_internal_error_ + call psb_errpush(info,name) + goto 9999 + end if + + ! + ! Compute destination for my data. + ! Simplistic reallocation: divide the NP processes + ! across the new ones (as balanced as possible), + ! then send all data from old to new process + ! + id2 = np/rnp + id1 = id2+1 + imd = mod(np,rnp) + if (me < (imd*id1)) then + ipdest = (me/id1) + else + ipdest = ( ((me-imd*id1)/id2) + imd) + end if + if (allocated(ids)) then + ipd = ids(ipdest) + else + ipd = ipdest + end if + write(0,*) ' Sending my data from ',me,' to ', & + & ipd, 'out of ',rnp,rnp-1 + + ! + ! Compute local rows for all new + ! processes; will have a BLOCK distribution + ! + newnl = 0 + newnl(ipdest+1) = desc_in%get_local_rows() + call psb_sum(ictxt,newnl) + + if (rme>=0) then + ! + if (rme < imd) then + isrc = [ (i, i=rme*id1,min(rme*id1+id1-1,np-1)) ] + else + isrc = [ (i, i= imd*id1+((rme-imd))*id2,& + & min(imd*id1+(rme-imd)*id2+id2-1,np-1) ) ] + end if + write(0,*) me,rme,imd,' ISRC: ',isrc(:) + nsrc = size(isrc) + write(0,*) me,rme,'In ',desc_in%get_local_rows(),desc_in%get_global_rows(),& + & ' out ',desc_out%get_local_rows(),desc_out%get_global_rows() + else + write(0,*) me,rme,'In ',desc_in%get_local_rows(),desc_in%get_global_rows(),& + & ' out ',0,0 + end if + else + write(0,*) 'Right now only BLOCK on input ' + info = psb_err_invalid_cd_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + ! + ! Collect matrices on their destinations + ! + block + integer(psb_ipk_) :: nzsnd, nzrcv, ip + integer(psb_ipk_) :: nrl, ncl, nzl, nzp + call a_in%cp_to(acoo_snd) + nzsnd = acoo_snd%get_nzeros() + call psb_snd(ictxt,nzsnd,ipd) + ! Convert to global numbering + call psb_loc_to_glob(acoo_snd%ia(1:nzsnd),desc_in,info) + call psb_loc_to_glob(acoo_snd%ja(1:nzsnd),desc_in,info) + + call psb_snd(ictxt,acoo_snd%ia(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%ja(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%val(1:nzsnd),ipd) + + if (rme>=0) then + ! prepare to receive + nzsrc = isrc + nzl = 0 + do ip=1, nsrc + call psb_rcv(ictxt,nzsrc(ip),isrc(ip)) + nzl = nzl + nzsrc(ip) + end do + call acoo_rcv%allocate(newnl(rme+1),newnl(rme+1),nzl) + nrl = acoo_rcv%get_nrows() + ncl = acoo_rcv%get_ncols() + nzp = 0 + do ip=1, nsrc + call psb_rcv(ictxt,acoo_rcv%ia(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + call psb_rcv(ictxt,acoo_rcv%ja(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + call psb_rcv(ictxt,acoo_rcv%val(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + nzp = nzp + nzsrc(ip) + end do + call acoo_rcv%set_nzeros(nzp) + write(0,*) rme,' Collected: ',& + & acoo_rcv%get_nrows(),acoo_rcv%get_ncols(),acoo_rcv%get_nzeros() + + ! + ! New descriptor + ! + call psb_cdall(newctxt,desc_out,info,nl=newnl(rme+1)) + ! Insert + call psb_spall(a_out,desc_out,info) + call psb_spins(nzp,acoo_rcv%ia(1:nzp),acoo_rcv%ja(1:nzp),& + & acoo_rcv%val(1:nzp),a_out,desc_out,info) + ! Assemble + call psb_cdasb(desc_out,info) + call psb_spasb(a_out,desc_out,info) + + write(0,*) rme,' Regenerated: ',& + & desc_out%get_local_rows(), desc_out%get_local_cols(),& + & a_out%get_nrows(),a_out%get_ncols(),a_out%get_nzeros() + +!!$ call desc_out%free(info) +!!$ call psb_exit(newctxt,close=.false.) + end if + end block + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_c_remap diff --git a/base/tools/psb_d_remap.F90 b/base/tools/psb_d_remap.F90 index 4835a0bc..e55626c5 100644 --- a/base/tools/psb_d_remap.F90 +++ b/base/tools/psb_d_remap.F90 @@ -30,8 +30,7 @@ ! ! ! -! Subroutine: psb_cdcpy -! Produces a clone of a descriptor. +! Subroutine: psb_d_remap ! ! Arguments: ! desc_in - type(psb_desc_type). The communication descriptor to be cloned. @@ -54,8 +53,8 @@ subroutine psb_d_remap(np_remap, desc_in, a_in, desc_out, a_out, info) ! locals integer(psb_ipk_) :: np, me, ictxt, err_act integer(psb_ipk_) :: rnp, rme, newctxt - integer(psb_ipk_) :: ipdest, id1, id2, imd, i, nsrc - integer(psb_ipk_), allocatable :: newnl(:), isrc(:), nzsrc(:) + integer(psb_ipk_) :: ipdest, ipd, id1, id2, imd, i, nsrc + integer(psb_ipk_), allocatable :: newnl(:), isrc(:), nzsrc(:), ids(:) type(psb_ld_coo_sparse_mat) :: acoo_snd, acoo_rcv integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -83,8 +82,28 @@ subroutine psb_d_remap(np_remap, desc_in, a_in, desc_out, a_out, info) write(0,*) ' Remapping from ',np,' onto ', np_remap if (desc_in%get_fmt() == 'BLOCK') then - ! OK - call psb_init(newctxt,np=np_remap,basectxt=ictxt) + ! + ! Should we spread the processes in the new context, + ! or should we keep them close? + ! + if (.true.) then + allocate(ids(0:np_remap-1)) + if (np_remap <= np/2) then + ids(0) = 0 + do ipdest=1,np_remap -1 + ids(ipdest) = ids(ipdest-1) + np/np_remap + end do + write(0,*) ' IDS ',ids(:) + else + do ipdest = 0, np_remap-1 + ids(ipdest) = ipdest + end do + end if + call psb_init(newctxt,np=np_remap,basectxt=ictxt,ids=ids) + else + call psb_init(newctxt,np=np_remap,basectxt=ictxt) + end if + call psb_info(newctxt,rme,rnp) write(0,*) 'Old context: ',me,np,' New context: ',rme,rnp call psb_bcast(ictxt,rnp) @@ -115,8 +134,13 @@ subroutine psb_d_remap(np_remap, desc_in, a_in, desc_out, a_out, info) else ipdest = ( ((me-imd*id1)/id2) + imd) end if + if (allocated(ids)) then + ipd = ids(ipdest) + else + ipd = ipdest + end if write(0,*) ' Sending my data from ',me,' to ', & - & ipdest, 'out of ',rnp,rnp-1 + & ipd, 'out of ',rnp,rnp-1 ! ! Compute local rows for all new @@ -158,14 +182,14 @@ subroutine psb_d_remap(np_remap, desc_in, a_in, desc_out, a_out, info) integer(psb_ipk_) :: nrl, ncl, nzl, nzp call a_in%cp_to(acoo_snd) nzsnd = acoo_snd%get_nzeros() - call psb_snd(ictxt,nzsnd,ipdest) + call psb_snd(ictxt,nzsnd,ipd) ! Convert to global numbering call psb_loc_to_glob(acoo_snd%ia(1:nzsnd),desc_in,info) call psb_loc_to_glob(acoo_snd%ja(1:nzsnd),desc_in,info) - call psb_snd(ictxt,acoo_snd%ia(1:nzsnd),ipdest) - call psb_snd(ictxt,acoo_snd%ja(1:nzsnd),ipdest) - call psb_snd(ictxt,acoo_snd%val(1:nzsnd),ipdest) + call psb_snd(ictxt,acoo_snd%ia(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%ja(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%val(1:nzsnd),ipd) if (rme>=0) then ! prepare to receive diff --git a/base/tools/psb_dasb.f90 b/base/tools/psb_dasb.f90 index 6aff978c..1d141cdf 100644 --- a/base/tools/psb_dasb.f90 +++ b/base/tools/psb_dasb.f90 @@ -80,25 +80,16 @@ subroutine psb_dasb_vect(x, desc_a, info, mold, scratch) scratch_ = .false. if (present(scratch)) scratch_ = scratch call psb_info(ictxt, me, np) - + ! ....verify blacs grid correctness.. - if (.false.) then - if (np == -1) then - info = psb_err_context_error_ - call psb_errpush(info,name) - goto 9999 - else if (.not.desc_a%is_ok()) then - info = psb_err_invalid_cd_state_ - call psb_errpush(info,name) - goto 9999 - end if - else - if (np == -1) then - ncol = 0 - call x%bld(ncol,mold=mold) - call psb_erractionrestore(err_act) - return - end if + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + else if (.not.desc_a%is_ok()) then + info = psb_err_invalid_cd_state_ + call psb_errpush(info,name) + goto 9999 end if nrow = desc_a%get_local_rows() diff --git a/base/tools/psb_s_remap.F90 b/base/tools/psb_s_remap.F90 new file mode 100644 index 00000000..6fa33223 --- /dev/null +++ b/base/tools/psb_s_remap.F90 @@ -0,0 +1,244 @@ +! +! 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. +! +! +! +! Subroutine: psb_s_remap +! +! Arguments: +! desc_in - type(psb_desc_type). The communication descriptor to be cloned. +! desc_out - type(psb_desc_type). The output communication descriptor. +! info - integer. Return code. +subroutine psb_s_remap(np_remap, desc_in, a_in, desc_out, a_out, info) + + use psb_base_mod, psb_protect_name => psb_s_remap + + implicit none + !....parameters... + integer(psb_ipk_), intent(in) :: np_remap + type(psb_desc_type), intent(inout) :: desc_in + type(psb_sspmat_type), intent(inout) :: a_in + type(psb_sspmat_type), intent(out) :: a_out + type(psb_desc_type), intent(out) :: desc_out + integer(psb_ipk_), intent(out) :: info + + + ! locals + integer(psb_ipk_) :: np, me, ictxt, err_act + integer(psb_ipk_) :: rnp, rme, newctxt + integer(psb_ipk_) :: ipdest, ipd, id1, id2, imd, i, nsrc + integer(psb_ipk_), allocatable :: newnl(:), isrc(:), nzsrc(:), ids(:) + type(psb_ls_coo_sparse_mat) :: acoo_snd, acoo_rcv + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name + + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if (psb_get_errstatus() /= 0) return + info = psb_success_ + call psb_erractionsave(err_act) + name = 'psb_cdcpy' + + ictxt = desc_in%get_context() + + ! check on blacs grid + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),': Entered' + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + write(0,*) ' Remapping from ',np,' onto ', np_remap + + if (desc_in%get_fmt() == 'BLOCK') then + ! + ! Should we spread the processes in the new context, + ! or should we keep them close? + ! + if (.true.) then + allocate(ids(0:np_remap-1)) + if (np_remap <= np/2) then + ids(0) = 0 + do ipdest=1,np_remap -1 + ids(ipdest) = ids(ipdest-1) + np/np_remap + end do + write(0,*) ' IDS ',ids(:) + else + do ipdest = 0, np_remap-1 + ids(ipdest) = ipdest + end do + end if + call psb_init(newctxt,np=np_remap,basectxt=ictxt,ids=ids) + else + call psb_init(newctxt,np=np_remap,basectxt=ictxt) + end if + + call psb_info(newctxt,rme,rnp) + write(0,*) 'Old context: ',me,np,' New context: ',rme,rnp + call psb_bcast(ictxt,rnp) + allocate(newnl(rnp),stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + if (rnp >= np) then + write(0,*) ' No remapping on larger proc count now' + info = psb_err_internal_error_ + call psb_errpush(info,name) + goto 9999 + end if + + ! + ! Compute destination for my data. + ! Simplistic reallocation: divide the NP processes + ! across the new ones (as balanced as possible), + ! then send all data from old to new process + ! + id2 = np/rnp + id1 = id2+1 + imd = mod(np,rnp) + if (me < (imd*id1)) then + ipdest = (me/id1) + else + ipdest = ( ((me-imd*id1)/id2) + imd) + end if + if (allocated(ids)) then + ipd = ids(ipdest) + else + ipd = ipdest + end if + write(0,*) ' Sending my data from ',me,' to ', & + & ipd, 'out of ',rnp,rnp-1 + + ! + ! Compute local rows for all new + ! processes; will have a BLOCK distribution + ! + newnl = 0 + newnl(ipdest+1) = desc_in%get_local_rows() + call psb_sum(ictxt,newnl) + + if (rme>=0) then + ! + if (rme < imd) then + isrc = [ (i, i=rme*id1,min(rme*id1+id1-1,np-1)) ] + else + isrc = [ (i, i= imd*id1+((rme-imd))*id2,& + & min(imd*id1+(rme-imd)*id2+id2-1,np-1) ) ] + end if + write(0,*) me,rme,imd,' ISRC: ',isrc(:) + nsrc = size(isrc) + write(0,*) me,rme,'In ',desc_in%get_local_rows(),desc_in%get_global_rows(),& + & ' out ',desc_out%get_local_rows(),desc_out%get_global_rows() + else + write(0,*) me,rme,'In ',desc_in%get_local_rows(),desc_in%get_global_rows(),& + & ' out ',0,0 + end if + else + write(0,*) 'Right now only BLOCK on input ' + info = psb_err_invalid_cd_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + ! + ! Collect matrices on their destinations + ! + block + integer(psb_ipk_) :: nzsnd, nzrcv, ip + integer(psb_ipk_) :: nrl, ncl, nzl, nzp + call a_in%cp_to(acoo_snd) + nzsnd = acoo_snd%get_nzeros() + call psb_snd(ictxt,nzsnd,ipd) + ! Convert to global numbering + call psb_loc_to_glob(acoo_snd%ia(1:nzsnd),desc_in,info) + call psb_loc_to_glob(acoo_snd%ja(1:nzsnd),desc_in,info) + + call psb_snd(ictxt,acoo_snd%ia(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%ja(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%val(1:nzsnd),ipd) + + if (rme>=0) then + ! prepare to receive + nzsrc = isrc + nzl = 0 + do ip=1, nsrc + call psb_rcv(ictxt,nzsrc(ip),isrc(ip)) + nzl = nzl + nzsrc(ip) + end do + call acoo_rcv%allocate(newnl(rme+1),newnl(rme+1),nzl) + nrl = acoo_rcv%get_nrows() + ncl = acoo_rcv%get_ncols() + nzp = 0 + do ip=1, nsrc + call psb_rcv(ictxt,acoo_rcv%ia(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + call psb_rcv(ictxt,acoo_rcv%ja(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + call psb_rcv(ictxt,acoo_rcv%val(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + nzp = nzp + nzsrc(ip) + end do + call acoo_rcv%set_nzeros(nzp) + write(0,*) rme,' Collected: ',& + & acoo_rcv%get_nrows(),acoo_rcv%get_ncols(),acoo_rcv%get_nzeros() + + ! + ! New descriptor + ! + call psb_cdall(newctxt,desc_out,info,nl=newnl(rme+1)) + ! Insert + call psb_spall(a_out,desc_out,info) + call psb_spins(nzp,acoo_rcv%ia(1:nzp),acoo_rcv%ja(1:nzp),& + & acoo_rcv%val(1:nzp),a_out,desc_out,info) + ! Assemble + call psb_cdasb(desc_out,info) + call psb_spasb(a_out,desc_out,info) + + write(0,*) rme,' Regenerated: ',& + & desc_out%get_local_rows(), desc_out%get_local_cols(),& + & a_out%get_nrows(),a_out%get_ncols(),a_out%get_nzeros() + +!!$ call desc_out%free(info) +!!$ call psb_exit(newctxt,close=.false.) + end if + end block + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_s_remap diff --git a/base/tools/psb_z_remap.F90 b/base/tools/psb_z_remap.F90 new file mode 100644 index 00000000..3802ec3d --- /dev/null +++ b/base/tools/psb_z_remap.F90 @@ -0,0 +1,244 @@ +! +! 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. +! +! +! +! Subroutine: psb_z_remap +! +! Arguments: +! desc_in - type(psb_desc_type). The communication descriptor to be cloned. +! desc_out - type(psb_desc_type). The output communication descriptor. +! info - integer. Return code. +subroutine psb_z_remap(np_remap, desc_in, a_in, desc_out, a_out, info) + + use psb_base_mod, psb_protect_name => psb_z_remap + + implicit none + !....parameters... + integer(psb_ipk_), intent(in) :: np_remap + type(psb_desc_type), intent(inout) :: desc_in + type(psb_zspmat_type), intent(inout) :: a_in + type(psb_zspmat_type), intent(out) :: a_out + type(psb_desc_type), intent(out) :: desc_out + integer(psb_ipk_), intent(out) :: info + + + ! locals + integer(psb_ipk_) :: np, me, ictxt, err_act + integer(psb_ipk_) :: rnp, rme, newctxt + integer(psb_ipk_) :: ipdest, ipd, id1, id2, imd, i, nsrc + integer(psb_ipk_), allocatable :: newnl(:), isrc(:), nzsrc(:), ids(:) + type(psb_lz_coo_sparse_mat) :: acoo_snd, acoo_rcv + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name + + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if (psb_get_errstatus() /= 0) return + info = psb_success_ + call psb_erractionsave(err_act) + name = 'psb_cdcpy' + + ictxt = desc_in%get_context() + + ! check on blacs grid + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_ext_) & + & write(debug_unit,*) me,' ',trim(name),': Entered' + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + write(0,*) ' Remapping from ',np,' onto ', np_remap + + if (desc_in%get_fmt() == 'BLOCK') then + ! + ! Should we spread the processes in the new context, + ! or should we keep them close? + ! + if (.true.) then + allocate(ids(0:np_remap-1)) + if (np_remap <= np/2) then + ids(0) = 0 + do ipdest=1,np_remap -1 + ids(ipdest) = ids(ipdest-1) + np/np_remap + end do + write(0,*) ' IDS ',ids(:) + else + do ipdest = 0, np_remap-1 + ids(ipdest) = ipdest + end do + end if + call psb_init(newctxt,np=np_remap,basectxt=ictxt,ids=ids) + else + call psb_init(newctxt,np=np_remap,basectxt=ictxt) + end if + + call psb_info(newctxt,rme,rnp) + write(0,*) 'Old context: ',me,np,' New context: ',rme,rnp + call psb_bcast(ictxt,rnp) + allocate(newnl(rnp),stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + if (rnp >= np) then + write(0,*) ' No remapping on larger proc count now' + info = psb_err_internal_error_ + call psb_errpush(info,name) + goto 9999 + end if + + ! + ! Compute destination for my data. + ! Simplistic reallocation: divide the NP processes + ! across the new ones (as balanced as possible), + ! then send all data from old to new process + ! + id2 = np/rnp + id1 = id2+1 + imd = mod(np,rnp) + if (me < (imd*id1)) then + ipdest = (me/id1) + else + ipdest = ( ((me-imd*id1)/id2) + imd) + end if + if (allocated(ids)) then + ipd = ids(ipdest) + else + ipd = ipdest + end if + write(0,*) ' Sending my data from ',me,' to ', & + & ipd, 'out of ',rnp,rnp-1 + + ! + ! Compute local rows for all new + ! processes; will have a BLOCK distribution + ! + newnl = 0 + newnl(ipdest+1) = desc_in%get_local_rows() + call psb_sum(ictxt,newnl) + + if (rme>=0) then + ! + if (rme < imd) then + isrc = [ (i, i=rme*id1,min(rme*id1+id1-1,np-1)) ] + else + isrc = [ (i, i= imd*id1+((rme-imd))*id2,& + & min(imd*id1+(rme-imd)*id2+id2-1,np-1) ) ] + end if + write(0,*) me,rme,imd,' ISRC: ',isrc(:) + nsrc = size(isrc) + write(0,*) me,rme,'In ',desc_in%get_local_rows(),desc_in%get_global_rows(),& + & ' out ',desc_out%get_local_rows(),desc_out%get_global_rows() + else + write(0,*) me,rme,'In ',desc_in%get_local_rows(),desc_in%get_global_rows(),& + & ' out ',0,0 + end if + else + write(0,*) 'Right now only BLOCK on input ' + info = psb_err_invalid_cd_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + ! + ! Collect matrices on their destinations + ! + block + integer(psb_ipk_) :: nzsnd, nzrcv, ip + integer(psb_ipk_) :: nrl, ncl, nzl, nzp + call a_in%cp_to(acoo_snd) + nzsnd = acoo_snd%get_nzeros() + call psb_snd(ictxt,nzsnd,ipd) + ! Convert to global numbering + call psb_loc_to_glob(acoo_snd%ia(1:nzsnd),desc_in,info) + call psb_loc_to_glob(acoo_snd%ja(1:nzsnd),desc_in,info) + + call psb_snd(ictxt,acoo_snd%ia(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%ja(1:nzsnd),ipd) + call psb_snd(ictxt,acoo_snd%val(1:nzsnd),ipd) + + if (rme>=0) then + ! prepare to receive + nzsrc = isrc + nzl = 0 + do ip=1, nsrc + call psb_rcv(ictxt,nzsrc(ip),isrc(ip)) + nzl = nzl + nzsrc(ip) + end do + call acoo_rcv%allocate(newnl(rme+1),newnl(rme+1),nzl) + nrl = acoo_rcv%get_nrows() + ncl = acoo_rcv%get_ncols() + nzp = 0 + do ip=1, nsrc + call psb_rcv(ictxt,acoo_rcv%ia(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + call psb_rcv(ictxt,acoo_rcv%ja(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + call psb_rcv(ictxt,acoo_rcv%val(nzp+1:nzp+nzsrc(ip)),isrc(ip)) + nzp = nzp + nzsrc(ip) + end do + call acoo_rcv%set_nzeros(nzp) + write(0,*) rme,' Collected: ',& + & acoo_rcv%get_nrows(),acoo_rcv%get_ncols(),acoo_rcv%get_nzeros() + + ! + ! New descriptor + ! + call psb_cdall(newctxt,desc_out,info,nl=newnl(rme+1)) + ! Insert + call psb_spall(a_out,desc_out,info) + call psb_spins(nzp,acoo_rcv%ia(1:nzp),acoo_rcv%ja(1:nzp),& + & acoo_rcv%val(1:nzp),a_out,desc_out,info) + ! Assemble + call psb_cdasb(desc_out,info) + call psb_spasb(a_out,desc_out,info) + + write(0,*) rme,' Regenerated: ',& + & desc_out%get_local_rows(), desc_out%get_local_cols(),& + & a_out%get_nrows(),a_out%get_ncols(),a_out%get_nzeros() + +!!$ call desc_out%free(info) +!!$ call psb_exit(newctxt,close=.false.) + end if + end block + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_z_remap