fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-06 22:55:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

Further changes for remote build, new method

Browse Source
This commit is contained in:
Salvatore Filippone
2022-02-03 13:25:51 +01:00
parent 0e676d2903
commit 6d0b26ecf1
30 changed files with 1498 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/internals/psi_graph_fnd_owner.F90
+7 -5
View File
@@ -41,6 +41,8 @@
!
! iprc(:) - integer(psb_ipk_), allocatable Output: process identifiers
! for the corresponding indices
! ladj(:) - integer(psb_ipk_), allocatable Output: A list of adjacent processes
!
! idxmap - class(psb_indx_map). The index map
! info - integer. return code.
!
@@ -76,7 +78,7 @@
! thereby limiting the memory footprint to a predefined maximum
! (that the user can force with psb_cd_set_maxspace()).
!
subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info)
subroutine psi_graph_fnd_owner(idx,iprc,ladj,idxmap,info)
use psb_serial_mod
use psb_const_mod
use psb_error_mod
@@ -93,13 +95,13 @@ subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info)
include 'mpif.h'
#endif
integer(psb_lpk_), intent(in) :: idx(:)
integer(psb_ipk_), allocatable, intent(out) :: iprc(:)
integer(psb_ipk_), allocatable, intent(out) :: iprc(:), ladj(:)
class(psb_indx_map), intent(inout) :: idxmap
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), allocatable :: tidx(:)
integer(psb_ipk_), allocatable :: tprc(:), tsmpl(:), ladj(:)
integer(psb_ipk_), allocatable :: tprc(:), tsmpl(:)
integer(psb_mpk_) :: icomm, minfo
integer(psb_ipk_) :: i,n_row,n_col,err_act,ip,j, nsampl_out,&
& nv, n_answers, nqries, nsampl_in, locr_max, ist, iend,&
@@ -208,7 +210,7 @@ subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info)
if (trace.and.(me == 0)) write(0,*) ' Initial sweep on user-defined topology',&
& nsampl_in
call psi_adj_fnd_sweep(idx,iprc,ladj,idxmap,nsampl_in,n_answers)
call idxmap%xtnd_p_adjcncy(ladj)
!call idxmap%xtnd_p_adjcncy(ladj)
nqries = nv - n_answers
nqries_max = nqries
call psb_max(ctxt,nqries_max)
@@ -259,7 +261,7 @@ subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info)
ladj = tprc(1:nlansw)
call psb_msort_unique(ladj,nadj)
call psb_realloc(nadj,ladj,info)
call idxmap%xtnd_p_adjcncy(ladj)
! call idxmap%xtnd_p_adjcncy(ladj)
if (do_timings) call psb_toc(idx_loop_a2a)
if (do_timings) call psb_tic(idx_loop_neigh)
!
base/internals/psi_indx_map_fnd_owner.F90
+5 -5
View File
@@ -72,7 +72,7 @@ subroutine psi_indx_map_fnd_owner(idx,iprc,idxmap,info)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), allocatable :: hhidx(:)
integer(psb_ipk_), allocatable :: hhidx(:), ladj(:)
integer(psb_mpk_) :: icomm, minfo
integer(psb_ipk_) :: i, err_act, hsize
integer(psb_lpk_) :: nv
@@ -183,7 +183,7 @@ subroutine psi_indx_map_fnd_owner(idx,iprc,idxmap,info)
tidx(k2) = idx(k1)
end if
end do
call psi_graph_fnd_owner(tidx,tprc,idxmap,info)
call psi_graph_fnd_owner(tidx,tprc,ladj,idxmap,info)
k2 = 0
do k1 = 1, nv
if (iprc(k1) < 0) then
@@ -198,10 +198,10 @@ subroutine psi_indx_map_fnd_owner(idx,iprc,idxmap,info)
end do
end block
else
call psi_graph_fnd_owner(idx,iprc,idxmap,info)
call psi_graph_fnd_owner(idx,iprc,ladj,idxmap,info)
end if
call idxmap%xtnd_p_adjcncy(ladj)
end if
if (gettime) then
base/modules/desc/psb_indx_map_mod.f90
+2 -1
View File
@@ -303,11 +303,12 @@ module psb_indx_map_mod
end interface
interface
subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info)
subroutine psi_graph_fnd_owner(idx,iprc,ladj,idxmap,info)
import :: psb_indx_map, psb_ipk_, psb_lpk_
implicit none
integer(psb_lpk_), intent(in) :: idx(:)
integer(psb_ipk_), allocatable, intent(out) :: iprc(:)
integer(psb_ipk_), allocatable, intent(out) :: ladj(:)
class(psb_indx_map), intent(inout) :: idxmap
integer(psb_ipk_), intent(out) :: info
end subroutine psi_graph_fnd_owner
base/modules/serial/psb_c_mat_mod.F90
+1 -1
View File
@@ -86,7 +86,7 @@ module psb_c_mat_mod
class(psb_c_base_sparse_mat), allocatable :: a
integer(psb_ipk_) :: remote_build=psb_matbld_noremote_
class(psb_lc_base_sparse_mat), allocatable :: rmta
type(psb_lc_coo_sparse_mat), allocatable :: rmta
contains
! Getters
base/modules/serial/psb_d_mat_mod.F90
+1 -1
View File
@@ -86,7 +86,7 @@ module psb_d_mat_mod
class(psb_d_base_sparse_mat), allocatable :: a
integer(psb_ipk_) :: remote_build=psb_matbld_noremote_
class(psb_ld_base_sparse_mat), allocatable :: rmta
type(psb_ld_coo_sparse_mat), allocatable :: rmta
contains
! Getters
base/modules/serial/psb_s_mat_mod.F90
+1 -1
View File
@@ -86,7 +86,7 @@ module psb_s_mat_mod
class(psb_s_base_sparse_mat), allocatable :: a
integer(psb_ipk_) :: remote_build=psb_matbld_noremote_
class(psb_ls_base_sparse_mat), allocatable :: rmta
type(psb_ls_coo_sparse_mat), allocatable :: rmta
contains
! Getters
base/modules/serial/psb_z_mat_mod.F90
+1 -1
View File
@@ -86,7 +86,7 @@ module psb_z_mat_mod
class(psb_z_base_sparse_mat), allocatable :: a
integer(psb_ipk_) :: remote_build=psb_matbld_noremote_
class(psb_lz_base_sparse_mat), allocatable :: rmta
type(psb_lz_coo_sparse_mat), allocatable :: rmta
contains
! Getters
base/modules/tools/psb_c_tools_mod.F90
+12 -1
View File
@@ -254,7 +254,7 @@ Module psb_c_tools_mod
import
implicit none
type(psb_cspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -262,6 +262,17 @@ Module psb_c_tools_mod
end subroutine psb_cspasb
end interface
interface psb_remote_mat
subroutine psb_lc_remote_mat(a,desc_a,b, info)
import
implicit none
type(psb_lc_coo_sparse_mat),Intent(inout) :: a
type(psb_desc_type),intent(inout) :: desc_a
type(psb_lc_coo_sparse_mat),Intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_lc_remote_mat
end interface psb_remote_mat
interface psb_spfree
subroutine psb_cspfree(a, desc_a,info)
import
base/modules/tools/psb_d_tools_mod.F90
+12 -1
View File
@@ -254,7 +254,7 @@ Module psb_d_tools_mod
import
implicit none
type(psb_dspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -262,6 +262,17 @@ Module psb_d_tools_mod
end subroutine psb_dspasb
end interface
interface psb_remote_mat
subroutine psb_ld_remote_mat(a,desc_a,b, info)
import
implicit none
type(psb_ld_coo_sparse_mat),Intent(inout) :: a
type(psb_desc_type),intent(inout) :: desc_a
type(psb_ld_coo_sparse_mat),Intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_ld_remote_mat
end interface psb_remote_mat
interface psb_spfree
subroutine psb_dspfree(a, desc_a,info)
import
base/modules/tools/psb_s_tools_mod.F90
+12 -1
View File
@@ -254,7 +254,7 @@ Module psb_s_tools_mod
import
implicit none
type(psb_sspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -262,6 +262,17 @@ Module psb_s_tools_mod
end subroutine psb_sspasb
end interface
interface psb_remote_mat
subroutine psb_ls_remote_mat(a,desc_a,b, info)
import
implicit none
type(psb_ls_coo_sparse_mat),Intent(inout) :: a
type(psb_desc_type),intent(inout) :: desc_a
type(psb_ls_coo_sparse_mat),Intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_ls_remote_mat
end interface psb_remote_mat
interface psb_spfree
subroutine psb_sspfree(a, desc_a,info)
import
base/modules/tools/psb_z_tools_mod.F90
+12 -1
View File
@@ -254,7 +254,7 @@ Module psb_z_tools_mod
import
implicit none
type(psb_zspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -262,6 +262,17 @@ Module psb_z_tools_mod
end subroutine psb_zspasb
end interface
interface psb_remote_mat
subroutine psb_lz_remote_mat(a,desc_a,b, info)
import
implicit none
type(psb_lz_coo_sparse_mat),Intent(inout) :: a
type(psb_desc_type),intent(inout) :: desc_a
type(psb_lz_coo_sparse_mat),Intent(inout) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_lz_remote_mat
end interface psb_remote_mat
interface psb_spfree
subroutine psb_zspfree(a, desc_a,info)
import
base/tools/Makefile
+2 -1
View File
@@ -30,7 +30,8 @@ FOBJS = psb_cdall.o psb_cdals.o psb_cdalv.o psb_cd_inloc.o psb_cdins.o psb_cdprt
psb_cgetelem.o psb_dgetelem.o psb_sgetelem.o psb_zgetelem.o
MPFOBJS = psb_icdasb.o psb_ssphalo.o psb_dsphalo.o psb_csphalo.o psb_zsphalo.o \
psb_dcdbldext.o psb_zcdbldext.o psb_scdbldext.o psb_ccdbldext.o
psb_dcdbldext.o psb_zcdbldext.o psb_scdbldext.o psb_ccdbldext.o \
psb_s_remote_mat.o psb_d_remote_mat.o psb_c_remote_mat.o psb_z_remote_mat.o
LIBDIR=..
INCDIR=..
base/tools/psb_c_remote_mat.F90
+259
View File
@@ -0,0 +1,259 @@
!
! 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: psb_csphalo.f90
!
! Subroutine: psb_csphalo psb_lcsphalo
! This routine does the retrieval of remote matrix rows.
! Retrieval is done through GETROW, therefore it should work
! for any matrix format in A; as for the output, default is CSR.
!
! There is also a specialized version lc_CSR whose interface
! is adapted for the needs of c_par_csr_spspmm.
!
! There are three possible exchange algorithms:
! 1. Use MPI_Alltoallv
! 2. Use psb_simple_a2av
! 3. Use psb_simple_triad_a2av
! Default choice is 3. The MPI variant has proved to be inefficient;
! that is because it is not persistent, therefore you pay the initialization price
! every time, and it is not optimized for a sparse communication pattern,
! most MPI implementations assume that all communications are non-empty.
! The PSB_SIMPLE variants reuse the same communicator, and go for a simplistic
! sequence of sends/receive that is quite efficient for a sparse communication
! pattern. To be refined/reviewed in the future to compare with neighbour
! persistent collectives.
!
! Arguments:
! a - type(psb_cspmat_type) The local part of input matrix A
! desc_a - type(psb_desc_type). The communication descriptor.
! blck - type(psb_cspmat_type) The local part of output matrix BLCK
! info - integer. Return code
! rowcnv - logical Should row/col indices be converted
! colcnv - logical to/from global numbering when sent/received?
! default is .TRUE.
! rowscale - logical Should row/col indices on output be remapped
! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ?
! default is .FALSE.
! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.)
! data - integer Which index list in desc_a should be used to retrieve
! rows, default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ DISABLED for this routine.
!
Subroutine psb_lc_remote_mat(a,desc_a,b,info)
use psb_base_mod, psb_protect_name => psb_lc_remote_mat
#ifdef MPI_MOD
use mpi
#endif
Implicit None
#ifdef MPI_H
include 'mpif.h'
#endif
Type(psb_lc_coo_sparse_mat),Intent(inout) :: a
Type(psb_lc_coo_sparse_mat),Intent(inout) :: b
Type(psb_desc_type), Intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
! ...local scalars....
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np,me
integer(psb_ipk_) :: counter, proc, i, n_el_send,n_el_recv, &
& n_elem, j, ipx,mat_recv, idxs,idxr,&
& data_,totxch,nxs, nxr, ncg
integer(psb_lpk_) :: r, k, irmin, irmax, icmin, icmax, iszs, iszr, &
& lidx, l1, lnr, lnc, lnnz, idx, ngtz, tot_elem
integer(psb_lpk_) :: nz,nouth
integer(psb_ipk_) :: nnp, nrcvs, nsnds
integer(psb_mpk_) :: icomm, minfo
integer(psb_mpk_), allocatable :: brvindx(:), &
& rvsz(:), bsdindx(:),sdsz(:)
integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:)
complex(psb_spk_), allocatable :: valsnd(:)
type(psb_lc_coo_sparse_mat), allocatable :: acoo
integer(psb_ipk_), pointer :: idxv(:)
class(psb_i_base_vect_type), pointer :: pdxv
integer(psb_ipk_), allocatable :: ipdxv(:), ladj(:), ila(:), iprc(:)
logical :: rowcnv_,colcnv_,rowscale_,colscale_
character(len=5) :: outfmt_
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
info=psb_success_
name='psb_csphalo'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
icomm = desc_a%get_mpic()
Call psb_info(ctxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
call b%free()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
nz = a%get_nzeros()
allocate(ila(nz))
write(0,*) me,name,' size :',nz,size(ila)
call desc_a%g2l(a%ia(1:nz),ila(1:nz),info,owned=.false.)
nouth = count(ila(1:nz)<0)
write(0,*) me,name,' Count out of halo :',nouth
call psb_max(ctxt,nouth)
if ((nouth/=0).and.(me==0)) &
& write(0,*) 'Warning: would require reinit of DESC_A'
call psi_graph_fnd_owner(a%ia(1:nz),iprc,ladj,desc_a%indxmap,info)
call psb_msort_unique(ladj,nnp)
write(0,*) me,name,' Processes:',ladj(1:nnp)
icomm = desc_a%get_mpic()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), stat=info)
sdsz(:)=0
rvsz(:)=0
ipx = 1
brvindx(ipx) = 0
bsdindx(ipx) = 0
counter=1
idx = 0
idxs = 0
idxr = 0
do i=1,nz
if (iprc(i) >=0) then
sdsz(iprc(i)+1) = sdsz(iprc(i)+1) +1
else
write(0,*)me,name,' Error from fnd_owner: ',iprc(i)
end if
end do
call mpi_alltoall(sdsz,1,psb_mpi_mpk_,&
& rvsz,1,psb_mpi_mpk_,icomm,minfo)
if (minfo /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='mpi_alltoall')
goto 9999
end if
write(0,*)me,name,' sdsz ',sdsz(:),' rvsz:',rvsz(:)
nsnds = count(sdsz /= 0)
nrcvs = count(rvsz /= 0)
idxs = 0
idxr = 0
counter = 1
lnnz = max(iszr,iszs,lone)
lnc = a%get_ncols()
call acoo%allocate(lnr,lnc,lnnz)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),&
& ' Send:',sdsz(:),' Receive:',rvsz(:)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='ensure_size')
goto 9999
end if
select case(psb_get_sp_a2av_alg())
case(psb_sp_a2av_smpl_triad_)
call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,&
& acoo%val,acoo%ia,acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_smpl_v_)
call psb_simple_a2av(valsnd,sdsz,bsdindx,&
& acoo%val,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,&
& acoo%ia,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,&
& acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_mpi_)
call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_c_spk_,&
& acoo%val,rvsz,brvindx,psb_mpi_c_spk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ia,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ja,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo /= mpi_success) info = minfo
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='wrong A2AV alg selector')
goto 9999
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='alltoallv')
goto 9999
end if
call acoo%mv_to_coo(b,info)
Deallocate(brvindx,bsdindx,rvsz,sdsz,&
& iasnd,jasnd,valsnd,stat=info)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Done'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
End Subroutine psb_lc_remote_mat
base/tools/psb_cspalloc.f90
+1 -1
View File
@@ -124,7 +124,7 @@ subroutine psb_cspalloc(a, desc_a, info, nnz, bldmode)
write(0,*) name,' matbld_noremote_ nothing needed'
case (psb_matbld_remote_)
write(0,*) name,' matbld_remote_ start '
allocate(psb_lc_coo_sparse_mat :: a%rmta)
allocate(a%rmta)
nnzrmt_ = max(100,(nnz_/100))
call a%rmta%allocate(m,n,nnzrmt_)
base/tools/psb_cspasb.f90
+10 -2
View File
@@ -50,13 +50,14 @@
!
subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold)
use psb_base_mod, psb_protect_name => psb_cspasb
use psb_sort_mod
use psi_mod
implicit none
!...Parameters....
type(psb_cspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -117,7 +118,14 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold)
case (psb_matbld_noremote_)
! nothing needed
case (psb_matbld_remote_)
write(0,*) me,' Size of rmta:',a%rmta%get_nzeros()
write(0,*) me,name,' Size of rmta:',a%rmta%get_nzeros()
block
type(psb_lc_coo_sparse_mat) :: a_add
call psb_remote_mat(a%rmta,desc_a,a_add,info)
end block
end select
call a%cscnv(info,type=afmt,dupl=dupl, mold=mold)
base/tools/psb_cspins.F90
+29 -3
View File
@@ -158,7 +158,7 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
write(0,*) 'Check on insert ',nnl
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
@@ -198,8 +198,9 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info)
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
& mask=(ila(1:nz)>0))
end if
call a%csput(nz,ila,jla,val,ione,nrow,ione,ncol,info)
if (info /= psb_success_) then
@@ -207,6 +208,31 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
call psb_errpush(info,name,a_err='a%csput')
goto 9999
end if
select case(a%remote_build)
case (psb_matbld_noremote_)
! Do nothing
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
if (ila(i)<0) then
k=k+1
lila(k) = ia(k)
ljla(k) = ja(k)
lval(k) = val(k)
end if
end do
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
case default
write(0,*) name,' Ignoring wrong value for %remote_build'
end select
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
base/tools/psb_d_remote_mat.F90
+259
View File
@@ -0,0 +1,259 @@
!
! 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: psb_dsphalo.f90
!
! Subroutine: psb_dsphalo psb_ldsphalo
! This routine does the retrieval of remote matrix rows.
! Retrieval is done through GETROW, therefore it should work
! for any matrix format in A; as for the output, default is CSR.
!
! There is also a specialized version ld_CSR whose interface
! is adapted for the needs of d_par_csr_spspmm.
!
! There are three possible exchange algorithms:
! 1. Use MPI_Alltoallv
! 2. Use psb_simple_a2av
! 3. Use psb_simple_triad_a2av
! Default choice is 3. The MPI variant has proved to be inefficient;
! that is because it is not persistent, therefore you pay the initialization price
! every time, and it is not optimized for a sparse communication pattern,
! most MPI implementations assume that all communications are non-empty.
! The PSB_SIMPLE variants reuse the same communicator, and go for a simplistic
! sequence of sends/receive that is quite efficient for a sparse communication
! pattern. To be refined/reviewed in the future to compare with neighbour
! persistent collectives.
!
! Arguments:
! a - type(psb_dspmat_type) The local part of input matrix A
! desc_a - type(psb_desc_type). The communication descriptor.
! blck - type(psb_dspmat_type) The local part of output matrix BLCK
! info - integer. Return code
! rowcnv - logical Should row/col indices be converted
! colcnv - logical to/from global numbering when sent/received?
! default is .TRUE.
! rowscale - logical Should row/col indices on output be remapped
! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ?
! default is .FALSE.
! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.)
! data - integer Which index list in desc_a should be used to retrieve
! rows, default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ DISABLED for this routine.
!
Subroutine psb_ld_remote_mat(a,desc_a,b,info)
use psb_base_mod, psb_protect_name => psb_ld_remote_mat
#ifdef MPI_MOD
use mpi
#endif
Implicit None
#ifdef MPI_H
include 'mpif.h'
#endif
Type(psb_ld_coo_sparse_mat),Intent(inout) :: a
Type(psb_ld_coo_sparse_mat),Intent(inout) :: b
Type(psb_desc_type), Intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
! ...local scalars....
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np,me
integer(psb_ipk_) :: counter, proc, i, n_el_send,n_el_recv, &
& n_elem, j, ipx,mat_recv, idxs,idxr,&
& data_,totxch,nxs, nxr, ncg
integer(psb_lpk_) :: r, k, irmin, irmax, icmin, icmax, iszs, iszr, &
& lidx, l1, lnr, lnc, lnnz, idx, ngtz, tot_elem
integer(psb_lpk_) :: nz,nouth
integer(psb_ipk_) :: nnp, nrcvs, nsnds
integer(psb_mpk_) :: icomm, minfo
integer(psb_mpk_), allocatable :: brvindx(:), &
& rvsz(:), bsdindx(:),sdsz(:)
integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:)
real(psb_dpk_), allocatable :: valsnd(:)
type(psb_ld_coo_sparse_mat), allocatable :: acoo
integer(psb_ipk_), pointer :: idxv(:)
class(psb_i_base_vect_type), pointer :: pdxv
integer(psb_ipk_), allocatable :: ipdxv(:), ladj(:), ila(:), iprc(:)
logical :: rowcnv_,colcnv_,rowscale_,colscale_
character(len=5) :: outfmt_
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
info=psb_success_
name='psb_dsphalo'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
icomm = desc_a%get_mpic()
Call psb_info(ctxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
call b%free()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
nz = a%get_nzeros()
allocate(ila(nz))
write(0,*) me,name,' size :',nz,size(ila)
call desc_a%g2l(a%ia(1:nz),ila(1:nz),info,owned=.false.)
nouth = count(ila(1:nz)<0)
write(0,*) me,name,' Count out of halo :',nouth
call psb_max(ctxt,nouth)
if ((nouth/=0).and.(me==0)) &
& write(0,*) 'Warning: would require reinit of DESC_A'
call psi_graph_fnd_owner(a%ia(1:nz),iprc,ladj,desc_a%indxmap,info)
call psb_msort_unique(ladj,nnp)
write(0,*) me,name,' Processes:',ladj(1:nnp)
icomm = desc_a%get_mpic()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), stat=info)
sdsz(:)=0
rvsz(:)=0
ipx = 1
brvindx(ipx) = 0
bsdindx(ipx) = 0
counter=1
idx = 0
idxs = 0
idxr = 0
do i=1,nz
if (iprc(i) >=0) then
sdsz(iprc(i)+1) = sdsz(iprc(i)+1) +1
else
write(0,*)me,name,' Error from fnd_owner: ',iprc(i)
end if
end do
call mpi_alltoall(sdsz,1,psb_mpi_mpk_,&
& rvsz,1,psb_mpi_mpk_,icomm,minfo)
if (minfo /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='mpi_alltoall')
goto 9999
end if
write(0,*)me,name,' sdsz ',sdsz(:),' rvsz:',rvsz(:)
nsnds = count(sdsz /= 0)
nrcvs = count(rvsz /= 0)
idxs = 0
idxr = 0
counter = 1
lnnz = max(iszr,iszs,lone)
lnc = a%get_ncols()
call acoo%allocate(lnr,lnc,lnnz)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),&
& ' Send:',sdsz(:),' Receive:',rvsz(:)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='ensure_size')
goto 9999
end if
select case(psb_get_sp_a2av_alg())
case(psb_sp_a2av_smpl_triad_)
call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,&
& acoo%val,acoo%ia,acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_smpl_v_)
call psb_simple_a2av(valsnd,sdsz,bsdindx,&
& acoo%val,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,&
& acoo%ia,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,&
& acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_mpi_)
call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_r_dpk_,&
& acoo%val,rvsz,brvindx,psb_mpi_r_dpk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ia,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ja,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo /= mpi_success) info = minfo
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='wrong A2AV alg selector')
goto 9999
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='alltoallv')
goto 9999
end if
call acoo%mv_to_coo(b,info)
Deallocate(brvindx,bsdindx,rvsz,sdsz,&
& iasnd,jasnd,valsnd,stat=info)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Done'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
End Subroutine psb_ld_remote_mat
base/tools/psb_dspalloc.f90
+1 -1
View File
@@ -124,7 +124,7 @@ subroutine psb_dspalloc(a, desc_a, info, nnz, bldmode)
write(0,*) name,' matbld_noremote_ nothing needed'
case (psb_matbld_remote_)
write(0,*) name,' matbld_remote_ start '
allocate(psb_ld_coo_sparse_mat :: a%rmta)
allocate(a%rmta)
nnzrmt_ = max(100,(nnz_/100))
call a%rmta%allocate(m,n,nnzrmt_)
base/tools/psb_dspasb.f90
+10 -2
View File
@@ -50,13 +50,14 @@
!
subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold)
use psb_base_mod, psb_protect_name => psb_dspasb
use psb_sort_mod
use psi_mod
implicit none
!...Parameters....
type(psb_dspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -117,7 +118,14 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold)
case (psb_matbld_noremote_)
! nothing needed
case (psb_matbld_remote_)
write(0,*) me,' Size of rmta:',a%rmta%get_nzeros()
write(0,*) me,name,' Size of rmta:',a%rmta%get_nzeros()
block
type(psb_ld_coo_sparse_mat) :: a_add
call psb_remote_mat(a%rmta,desc_a,a_add,info)
end block
end select
call a%cscnv(info,type=afmt,dupl=dupl, mold=mold)
base/tools/psb_dspins.F90
+29 -4
View File
@@ -156,10 +156,9 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
case (psb_matbld_noremote_)
! Do nothing
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
write(0,*) 'Check on insert ',nnl
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
@@ -199,8 +198,9 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info)
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
& mask=(ila(1:nz)>0))
end if
call a%csput(nz,ila,jla,val,ione,nrow,ione,ncol,info)
if (info /= psb_success_) then
@@ -208,6 +208,31 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
call psb_errpush(info,name,a_err='a%csput')
goto 9999
end if
select case(a%remote_build)
case (psb_matbld_noremote_)
! Do nothing
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
if (ila(i)<0) then
k=k+1
lila(k) = ia(k)
ljla(k) = ja(k)
lval(k) = val(k)
end if
end do
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
case default
write(0,*) name,' Ignoring wrong value for %remote_build'
end select
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
base/tools/psb_s_remote_mat.F90
+259
View File
@@ -0,0 +1,259 @@
!
! 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: psb_ssphalo.f90
!
! Subroutine: psb_ssphalo psb_lssphalo
! This routine does the retrieval of remote matrix rows.
! Retrieval is done through GETROW, therefore it should work
! for any matrix format in A; as for the output, default is CSR.
!
! There is also a specialized version ls_CSR whose interface
! is adapted for the needs of s_par_csr_spspmm.
!
! There are three possible exchange algorithms:
! 1. Use MPI_Alltoallv
! 2. Use psb_simple_a2av
! 3. Use psb_simple_triad_a2av
! Default choice is 3. The MPI variant has proved to be inefficient;
! that is because it is not persistent, therefore you pay the initialization price
! every time, and it is not optimized for a sparse communication pattern,
! most MPI implementations assume that all communications are non-empty.
! The PSB_SIMPLE variants reuse the same communicator, and go for a simplistic
! sequence of sends/receive that is quite efficient for a sparse communication
! pattern. To be refined/reviewed in the future to compare with neighbour
! persistent collectives.
!
! Arguments:
! a - type(psb_sspmat_type) The local part of input matrix A
! desc_a - type(psb_desc_type). The communication descriptor.
! blck - type(psb_sspmat_type) The local part of output matrix BLCK
! info - integer. Return code
! rowcnv - logical Should row/col indices be converted
! colcnv - logical to/from global numbering when sent/received?
! default is .TRUE.
! rowscale - logical Should row/col indices on output be remapped
! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ?
! default is .FALSE.
! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.)
! data - integer Which index list in desc_a should be used to retrieve
! rows, default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ DISABLED for this routine.
!
Subroutine psb_ls_remote_mat(a,desc_a,b,info)
use psb_base_mod, psb_protect_name => psb_ls_remote_mat
#ifdef MPI_MOD
use mpi
#endif
Implicit None
#ifdef MPI_H
include 'mpif.h'
#endif
Type(psb_ls_coo_sparse_mat),Intent(inout) :: a
Type(psb_ls_coo_sparse_mat),Intent(inout) :: b
Type(psb_desc_type), Intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
! ...local scalars....
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np,me
integer(psb_ipk_) :: counter, proc, i, n_el_send,n_el_recv, &
& n_elem, j, ipx,mat_recv, idxs,idxr,&
& data_,totxch,nxs, nxr, ncg
integer(psb_lpk_) :: r, k, irmin, irmax, icmin, icmax, iszs, iszr, &
& lidx, l1, lnr, lnc, lnnz, idx, ngtz, tot_elem
integer(psb_lpk_) :: nz,nouth
integer(psb_ipk_) :: nnp, nrcvs, nsnds
integer(psb_mpk_) :: icomm, minfo
integer(psb_mpk_), allocatable :: brvindx(:), &
& rvsz(:), bsdindx(:),sdsz(:)
integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:)
real(psb_spk_), allocatable :: valsnd(:)
type(psb_ls_coo_sparse_mat), allocatable :: acoo
integer(psb_ipk_), pointer :: idxv(:)
class(psb_i_base_vect_type), pointer :: pdxv
integer(psb_ipk_), allocatable :: ipdxv(:), ladj(:), ila(:), iprc(:)
logical :: rowcnv_,colcnv_,rowscale_,colscale_
character(len=5) :: outfmt_
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
info=psb_success_
name='psb_ssphalo'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
icomm = desc_a%get_mpic()
Call psb_info(ctxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
call b%free()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
nz = a%get_nzeros()
allocate(ila(nz))
write(0,*) me,name,' size :',nz,size(ila)
call desc_a%g2l(a%ia(1:nz),ila(1:nz),info,owned=.false.)
nouth = count(ila(1:nz)<0)
write(0,*) me,name,' Count out of halo :',nouth
call psb_max(ctxt,nouth)
if ((nouth/=0).and.(me==0)) &
& write(0,*) 'Warning: would require reinit of DESC_A'
call psi_graph_fnd_owner(a%ia(1:nz),iprc,ladj,desc_a%indxmap,info)
call psb_msort_unique(ladj,nnp)
write(0,*) me,name,' Processes:',ladj(1:nnp)
icomm = desc_a%get_mpic()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), stat=info)
sdsz(:)=0
rvsz(:)=0
ipx = 1
brvindx(ipx) = 0
bsdindx(ipx) = 0
counter=1
idx = 0
idxs = 0
idxr = 0
do i=1,nz
if (iprc(i) >=0) then
sdsz(iprc(i)+1) = sdsz(iprc(i)+1) +1
else
write(0,*)me,name,' Error from fnd_owner: ',iprc(i)
end if
end do
call mpi_alltoall(sdsz,1,psb_mpi_mpk_,&
& rvsz,1,psb_mpi_mpk_,icomm,minfo)
if (minfo /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='mpi_alltoall')
goto 9999
end if
write(0,*)me,name,' sdsz ',sdsz(:),' rvsz:',rvsz(:)
nsnds = count(sdsz /= 0)
nrcvs = count(rvsz /= 0)
idxs = 0
idxr = 0
counter = 1
lnnz = max(iszr,iszs,lone)
lnc = a%get_ncols()
call acoo%allocate(lnr,lnc,lnnz)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),&
& ' Send:',sdsz(:),' Receive:',rvsz(:)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='ensure_size')
goto 9999
end if
select case(psb_get_sp_a2av_alg())
case(psb_sp_a2av_smpl_triad_)
call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,&
& acoo%val,acoo%ia,acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_smpl_v_)
call psb_simple_a2av(valsnd,sdsz,bsdindx,&
& acoo%val,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,&
& acoo%ia,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,&
& acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_mpi_)
call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_r_spk_,&
& acoo%val,rvsz,brvindx,psb_mpi_r_spk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ia,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ja,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo /= mpi_success) info = minfo
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='wrong A2AV alg selector')
goto 9999
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='alltoallv')
goto 9999
end if
call acoo%mv_to_coo(b,info)
Deallocate(brvindx,bsdindx,rvsz,sdsz,&
& iasnd,jasnd,valsnd,stat=info)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Done'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
End Subroutine psb_ls_remote_mat
base/tools/psb_sspalloc.f90
+1 -1
View File
@@ -124,7 +124,7 @@ subroutine psb_sspalloc(a, desc_a, info, nnz, bldmode)
write(0,*) name,' matbld_noremote_ nothing needed'
case (psb_matbld_remote_)
write(0,*) name,' matbld_remote_ start '
allocate(psb_ls_coo_sparse_mat :: a%rmta)
allocate(a%rmta)
nnzrmt_ = max(100,(nnz_/100))
call a%rmta%allocate(m,n,nnzrmt_)
base/tools/psb_sspasb.f90
+10 -2
View File
@@ -50,13 +50,14 @@
!
subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold)
use psb_base_mod, psb_protect_name => psb_sspasb
use psb_sort_mod
use psi_mod
implicit none
!...Parameters....
type(psb_sspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -117,7 +118,14 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold)
case (psb_matbld_noremote_)
! nothing needed
case (psb_matbld_remote_)
write(0,*) me,' Size of rmta:',a%rmta%get_nzeros()
write(0,*) me,name,' Size of rmta:',a%rmta%get_nzeros()
block
type(psb_ls_coo_sparse_mat) :: a_add
call psb_remote_mat(a%rmta,desc_a,a_add,info)
end block
end select
call a%cscnv(info,type=afmt,dupl=dupl, mold=mold)
base/tools/psb_sspins.F90
+29 -3
View File
@@ -158,7 +158,7 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
write(0,*) 'Check on insert ',nnl
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
@@ -198,8 +198,9 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info)
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
& mask=(ila(1:nz)>0))
end if
call a%csput(nz,ila,jla,val,ione,nrow,ione,ncol,info)
if (info /= psb_success_) then
@@ -207,6 +208,31 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
call psb_errpush(info,name,a_err='a%csput')
goto 9999
end if
select case(a%remote_build)
case (psb_matbld_noremote_)
! Do nothing
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
if (ila(i)<0) then
k=k+1
lila(k) = ia(k)
ljla(k) = ja(k)
lval(k) = val(k)
end if
end do
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
case default
write(0,*) name,' Ignoring wrong value for %remote_build'
end select
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
base/tools/psb_z_remote_mat.F90
+259
View File
@@ -0,0 +1,259 @@
!
! 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: psb_zsphalo.f90
!
! Subroutine: psb_zsphalo psb_lzsphalo
! This routine does the retrieval of remote matrix rows.
! Retrieval is done through GETROW, therefore it should work
! for any matrix format in A; as for the output, default is CSR.
!
! There is also a specialized version lz_CSR whose interface
! is adapted for the needs of z_par_csr_spspmm.
!
! There are three possible exchange algorithms:
! 1. Use MPI_Alltoallv
! 2. Use psb_simple_a2av
! 3. Use psb_simple_triad_a2av
! Default choice is 3. The MPI variant has proved to be inefficient;
! that is because it is not persistent, therefore you pay the initialization price
! every time, and it is not optimized for a sparse communication pattern,
! most MPI implementations assume that all communications are non-empty.
! The PSB_SIMPLE variants reuse the same communicator, and go for a simplistic
! sequence of sends/receive that is quite efficient for a sparse communication
! pattern. To be refined/reviewed in the future to compare with neighbour
! persistent collectives.
!
! Arguments:
! a - type(psb_zspmat_type) The local part of input matrix A
! desc_a - type(psb_desc_type). The communication descriptor.
! blck - type(psb_zspmat_type) The local part of output matrix BLCK
! info - integer. Return code
! rowcnv - logical Should row/col indices be converted
! colcnv - logical to/from global numbering when sent/received?
! default is .TRUE.
! rowscale - logical Should row/col indices on output be remapped
! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ?
! default is .FALSE.
! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.)
! data - integer Which index list in desc_a should be used to retrieve
! rows, default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ DISABLED for this routine.
!
Subroutine psb_lz_remote_mat(a,desc_a,b,info)
use psb_base_mod, psb_protect_name => psb_lz_remote_mat
#ifdef MPI_MOD
use mpi
#endif
Implicit None
#ifdef MPI_H
include 'mpif.h'
#endif
Type(psb_lz_coo_sparse_mat),Intent(inout) :: a
Type(psb_lz_coo_sparse_mat),Intent(inout) :: b
Type(psb_desc_type), Intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
! ...local scalars....
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np,me
integer(psb_ipk_) :: counter, proc, i, n_el_send,n_el_recv, &
& n_elem, j, ipx,mat_recv, idxs,idxr,&
& data_,totxch,nxs, nxr, ncg
integer(psb_lpk_) :: r, k, irmin, irmax, icmin, icmax, iszs, iszr, &
& lidx, l1, lnr, lnc, lnnz, idx, ngtz, tot_elem
integer(psb_lpk_) :: nz,nouth
integer(psb_ipk_) :: nnp, nrcvs, nsnds
integer(psb_mpk_) :: icomm, minfo
integer(psb_mpk_), allocatable :: brvindx(:), &
& rvsz(:), bsdindx(:),sdsz(:)
integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:)
complex(psb_dpk_), allocatable :: valsnd(:)
type(psb_lz_coo_sparse_mat), allocatable :: acoo
integer(psb_ipk_), pointer :: idxv(:)
class(psb_i_base_vect_type), pointer :: pdxv
integer(psb_ipk_), allocatable :: ipdxv(:), ladj(:), ila(:), iprc(:)
logical :: rowcnv_,colcnv_,rowscale_,colscale_
character(len=5) :: outfmt_
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
info=psb_success_
name='psb_zsphalo'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
icomm = desc_a%get_mpic()
Call psb_info(ctxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
call b%free()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
nz = a%get_nzeros()
allocate(ila(nz))
write(0,*) me,name,' size :',nz,size(ila)
call desc_a%g2l(a%ia(1:nz),ila(1:nz),info,owned=.false.)
nouth = count(ila(1:nz)<0)
write(0,*) me,name,' Count out of halo :',nouth
call psb_max(ctxt,nouth)
if ((nouth/=0).and.(me==0)) &
& write(0,*) 'Warning: would require reinit of DESC_A'
call psi_graph_fnd_owner(a%ia(1:nz),iprc,ladj,desc_a%indxmap,info)
call psb_msort_unique(ladj,nnp)
write(0,*) me,name,' Processes:',ladj(1:nnp)
icomm = desc_a%get_mpic()
Allocate(brvindx(np+1),&
& rvsz(np),sdsz(np),bsdindx(np+1), stat=info)
sdsz(:)=0
rvsz(:)=0
ipx = 1
brvindx(ipx) = 0
bsdindx(ipx) = 0
counter=1
idx = 0
idxs = 0
idxr = 0
do i=1,nz
if (iprc(i) >=0) then
sdsz(iprc(i)+1) = sdsz(iprc(i)+1) +1
else
write(0,*)me,name,' Error from fnd_owner: ',iprc(i)
end if
end do
call mpi_alltoall(sdsz,1,psb_mpi_mpk_,&
& rvsz,1,psb_mpi_mpk_,icomm,minfo)
if (minfo /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='mpi_alltoall')
goto 9999
end if
write(0,*)me,name,' sdsz ',sdsz(:),' rvsz:',rvsz(:)
nsnds = count(sdsz /= 0)
nrcvs = count(rvsz /= 0)
idxs = 0
idxr = 0
counter = 1
lnnz = max(iszr,iszs,lone)
lnc = a%get_ncols()
call acoo%allocate(lnr,lnc,lnnz)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),&
& ' Send:',sdsz(:),' Receive:',rvsz(:)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info)
if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='ensure_size')
goto 9999
end if
select case(psb_get_sp_a2av_alg())
case(psb_sp_a2av_smpl_triad_)
call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,&
& acoo%val,acoo%ia,acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_smpl_v_)
call psb_simple_a2av(valsnd,sdsz,bsdindx,&
& acoo%val,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,&
& acoo%ia,rvsz,brvindx,ctxt,info)
if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,&
& acoo%ja,rvsz,brvindx,ctxt,info)
case(psb_sp_a2av_mpi_)
call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_c_dpk_,&
& acoo%val,rvsz,brvindx,psb_mpi_c_dpk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ia,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo == mpi_success) &
& call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,&
& acoo%ja,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
if (minfo /= mpi_success) info = minfo
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='wrong A2AV alg selector')
goto 9999
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='alltoallv')
goto 9999
end if
call acoo%mv_to_coo(b,info)
Deallocate(brvindx,bsdindx,rvsz,sdsz,&
& iasnd,jasnd,valsnd,stat=info)
if (debug_level >= psb_debug_outer_)&
& write(debug_unit,*) me,' ',trim(name),': Done'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ctxt,err_act)
return
End Subroutine psb_lz_remote_mat
base/tools/psb_zspalloc.f90
+1 -1
View File
@@ -124,7 +124,7 @@ subroutine psb_zspalloc(a, desc_a, info, nnz, bldmode)
write(0,*) name,' matbld_noremote_ nothing needed'
case (psb_matbld_remote_)
write(0,*) name,' matbld_remote_ start '
allocate(psb_lz_coo_sparse_mat :: a%rmta)
allocate(a%rmta)
nnzrmt_ = max(100,(nnz_/100))
call a%rmta%allocate(m,n,nnzrmt_)
base/tools/psb_zspasb.f90
+10 -2
View File
@@ -50,13 +50,14 @@
!
subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold)
use psb_base_mod, psb_protect_name => psb_zspasb
use psb_sort_mod
use psi_mod
implicit none
!...Parameters....
type(psb_zspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(inout) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
character(len=*), optional, intent(in) :: afmt
@@ -117,7 +118,14 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold)
case (psb_matbld_noremote_)
! nothing needed
case (psb_matbld_remote_)
write(0,*) me,' Size of rmta:',a%rmta%get_nzeros()
write(0,*) me,name,' Size of rmta:',a%rmta%get_nzeros()
block
type(psb_lz_coo_sparse_mat) :: a_add
call psb_remote_mat(a%rmta,desc_a,a_add,info)
end block
end select
call a%cscnv(info,type=afmt,dupl=dupl, mold=mold)
base/tools/psb_zspins.F90
+29 -3
View File
@@ -158,7 +158,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
write(0,*) 'Check on insert ',nnl
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
@@ -198,8 +198,9 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info)
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
& mask=(ila(1:nz)>0))
end if
call a%csput(nz,ila,jla,val,ione,nrow,ione,ncol,info)
if (info /= psb_success_) then
@@ -207,6 +208,31 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
call psb_errpush(info,name,a_err='a%csput')
goto 9999
end if
select case(a%remote_build)
case (psb_matbld_noremote_)
! Do nothing
case (psb_matbld_remote_)
nnl = count(ila(1:nz)<0)
if (nnl > 0) then
!write(0,*) 'Check on insert ',nnl
allocate(lila(nnl),ljla(nnl),lval(nnl))
k = 0
do i=1,nz
if (ila(i)<0) then
k=k+1
lila(k) = ia(k)
ljla(k) = ja(k)
lval(k) = val(k)
end if
end do
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
case default
write(0,*) name,' Ignoring wrong value for %remote_build'
end select
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)
test/pargen/psb_d_pde3d.F90
+113 -28
View File
@@ -377,30 +377,7 @@ contains
!
call psb_cdall(ctxt,desc_a,info,vl=myidx)
!
! Add extra rows
!
block
integer(psb_ipk_) :: ks
mysz = nlr
if (m>nlr) mysz = mysz + m/nlr
call psb_realloc(mysz,myidx,info)
ks = nlr
outer: do i=1,idim
do j=1,idim
do k=1,idim
if (outside(i,j,k,bndx,bndy,bndz,iamx,iamy,iamz)) then
ks = ks + 1
if (ks > mysz) exit outer
call ijk2idx(myidx(ks),i,j,k,idim,idim,idim)
end if
end do
end do
end do outer
write(0,*) iam,' Check on extra nodes ',nlr,mysz,':',myidx(nlr+1:mysz)
end block
!
! Specify process topology
!
@@ -487,9 +464,9 @@ contains
call psb_barrier(ctxt)
t1 = psb_wtime()
do ii=1, mysz, nb
!ib = min(nb,nlr-ii+1)
ib = min(nb,mysz-ii+1)
do ii=1, nlr, nb
ib = min(nb,nlr-ii+1)
!ib = min(nb,mysz-ii+1)
icoeff = 1
do k=1,ib
i=ii+k-1
@@ -585,12 +562,119 @@ contains
goto 9999
end if
deallocate(val,irow,icol)
call psb_barrier(ctxt)
t1 = psb_wtime()
call psb_cdasb(desc_a,info,mold=imold)
tcdasb = psb_wtime()-t1
!
! Add extra rows
!
block
integer(psb_ipk_) :: ks, i
ks = desc_a%get_local_cols()-desc_a%get_local_rows()
if (ks > 0) ks = max(1,ks / 10)
mysz = nlr+ks
call psb_realloc(mysz,myidx,info)
do i=nlr+1, mysz
myidx(i) = i
end do
call desc_a%l2gv1(myidx(nlr+1:mysz),info)
!write(0,*) iam,' Check on extra nodes ',nlr,mysz,':',myidx(nlr+1:mysz)
do ii= nlr+1, mysz, nb
ib = min(nb,mysz-ii+1)
icoeff = 1
do k=1,ib
i=ii+k-1
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
zt(k) = f_(x,y,z)
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
if (ix == 1) then
zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y-1,z)
val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
if (iy == 1) then
zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z-1)
val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
if (iz == 1) then
zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z)
val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
& + c(x,y,z)
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y,z+1)
val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
if (iz == idim) then
zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y+1,z)
val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
if (iy == idim) then
zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x+1,y,z)
val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
if (ix==idim) then
zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
end do
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=dzero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
end block
call psb_barrier(ctxt)
t1 = psb_wtime()
if (info == psb_success_) then
@@ -634,6 +718,7 @@ contains
write(psb_out_unit,'("-total time : ",es12.5)') ttot
end if
deallocate(val,irow,icol)
call psb_erractionrestore(err_act)
return
test/pargen/psb_s_pde3d.F90
+121 -4
View File
@@ -208,7 +208,7 @@ contains
type(psb_s_coo_sparse_mat) :: acoo
type(psb_s_csr_sparse_mat) :: acsr
real(psb_spk_) :: zt(nb),x,y,z
integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_, mysz
integer(psb_lpk_) :: m,n,glob_row,nt
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
! For 3D partition
@@ -377,6 +377,7 @@ contains
!
call psb_cdall(ctxt,desc_a,info,vl=myidx)
!
! Specify process topology
!
@@ -463,8 +464,9 @@ contains
call psb_barrier(ctxt)
t1 = psb_wtime()
do ii=1, nlr,nb
do ii=1, nlr, nb
ib = min(nb,nlr-ii+1)
!ib = min(nb,mysz-ii+1)
icoeff = 1
do k=1,ib
i=ii+k-1
@@ -560,12 +562,119 @@ contains
goto 9999
end if
deallocate(val,irow,icol)
call psb_barrier(ctxt)
t1 = psb_wtime()
call psb_cdasb(desc_a,info,mold=imold)
tcdasb = psb_wtime()-t1
!
! Add extra rows
!
block
integer(psb_ipk_) :: ks, i
ks = desc_a%get_local_cols()-desc_a%get_local_rows()
if (ks > 0) ks = max(1,ks / 10)
mysz = nlr+ks
call psb_realloc(mysz,myidx,info)
do i=nlr+1, mysz
myidx(i) = i
end do
call desc_a%l2gv1(myidx(nlr+1:mysz),info)
!write(0,*) iam,' Check on extra nodes ',nlr,mysz,':',myidx(nlr+1:mysz)
do ii= nlr+1, mysz, nb
ib = min(nb,mysz-ii+1)
icoeff = 1
do k=1,ib
i=ii+k-1
! local matrix pointer
glob_row=myidx(i)
! compute gridpoint coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
zt(k) = f_(x,y,z)
! internal point: build discretization
!
! term depending on (x-1,y,z)
!
val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
if (ix == 1) then
zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y-1,z)
val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
if (iy == 1) then
zt(k) = g(x,szero,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z-1)
val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
if (iz == 1) then
zt(k) = g(x,y,szero)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y,z)
val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
& + c(x,y,z)
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
! term depending on (x,y,z+1)
val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
if (iz == idim) then
zt(k) = g(x,y,sone)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x,y+1,z)
val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
if (iy == idim) then
zt(k) = g(x,sone,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
! term depending on (x+1,y,z)
val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
if (ix==idim) then
zt(k) = g(sone,y,z)*(-val(icoeff)) + zt(k)
else
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
irow(icoeff) = glob_row
icoeff = icoeff+1
endif
end do
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
if(info /= psb_success_) exit
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
if(info /= psb_success_) exit
zt(:)=szero
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
if(info /= psb_success_) exit
end do
end block
call psb_barrier(ctxt)
t1 = psb_wtime()
if (info == psb_success_) then
@@ -609,6 +718,7 @@ contains
write(psb_out_unit,'("-total time : ",es12.5)') ttot
end if
deallocate(val,irow,icol)
call psb_erractionrestore(err_act)
return
@@ -616,8 +726,15 @@ contains
return
end subroutine psb_s_gen_pde3d
function outside(i,j,k,bndx,bndy,bndz,iamx,iamy,iamz) result(res)
logical :: res
integer(psb_ipk_), intent(in) :: i,j,k,iamx,iamy,iamz
integer(psb_ipk_), intent(in) :: bndx(0:),bndy(0:),bndz(0:)
res = (i<bndx(iamx)).or.(i>=bndx(iamx+1)) &
& .or.(j<bndy(iamy)).or.(j>=bndy(iamy+1)) &
& .or.(k<bndz(iamz)).or.(k>=bndz(iamz+1))
end function outside
end module psb_s_pde3d_mod
program psb_s_pde3d