From 1337009f910ba48f0a84f24144bc558101ad1be1 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 1 Feb 2022 10:56:06 +0100 Subject: [PATCH 01/28] Infrastructure for non-local matrix build --- base/modules/desc/psb_hash_map_mod.f90 | 3 ++- base/modules/psb_const_mod.F90 | 3 +++ base/modules/serial/psb_c_mat_mod.F90 | 2 ++ base/modules/serial/psb_d_mat_mod.F90 | 2 ++ base/modules/serial/psb_s_mat_mod.F90 | 2 ++ base/modules/serial/psb_z_mat_mod.F90 | 2 ++ base/modules/tools/psb_c_tools_mod.F90 | 4 +-- base/modules/tools/psb_d_tools_mod.F90 | 4 +-- base/modules/tools/psb_s_tools_mod.F90 | 4 +-- base/modules/tools/psb_z_tools_mod.F90 | 4 +-- base/serial/impl/psb_c_coo_impl.F90 | 1 - base/serial/impl/psb_c_mat_impl.F90 | 7 +++++- base/serial/impl/psb_d_coo_impl.F90 | 1 - base/serial/impl/psb_d_mat_impl.F90 | 7 +++++- base/serial/impl/psb_s_coo_impl.F90 | 1 - base/serial/impl/psb_s_mat_impl.F90 | 7 +++++- base/serial/impl/psb_z_coo_impl.F90 | 1 - base/serial/impl/psb_z_mat_impl.F90 | 7 +++++- base/tools/psb_cspalloc.f90 | 34 ++++++++++++++++++++++---- base/tools/psb_cspasb.f90 | 10 +++++++- base/tools/psb_cspins.F90 | 28 +++++++++++++++++++++ base/tools/psb_dspalloc.f90 | 34 ++++++++++++++++++++++---- base/tools/psb_dspasb.f90 | 10 +++++++- base/tools/psb_dspins.F90 | 28 +++++++++++++++++++++ base/tools/psb_sspalloc.f90 | 34 ++++++++++++++++++++++---- base/tools/psb_sspasb.f90 | 10 +++++++- base/tools/psb_sspins.F90 | 28 +++++++++++++++++++++ base/tools/psb_zspalloc.f90 | 34 ++++++++++++++++++++++---- base/tools/psb_zspasb.f90 | 10 +++++++- base/tools/psb_zspins.F90 | 28 +++++++++++++++++++++ cbind/base/psb_c_serial_cbind_mod.F90 | 33 ------------------------- test/pargen/psb_d_pde3d.F90 | 3 ++- test/pargen/psb_s_pde3d.F90 | 3 ++- 33 files changed, 313 insertions(+), 76 deletions(-) diff --git a/base/modules/desc/psb_hash_map_mod.f90 b/base/modules/desc/psb_hash_map_mod.f90 index 3cfd33a4..5fe7b75b 100644 --- a/base/modules/desc/psb_hash_map_mod.f90 +++ b/base/modules/desc/psb_hash_map_mod.f90 @@ -746,7 +746,8 @@ contains if (info >=0) then if (nxt == lip) then ncol = max(nxt,ncol) - call psb_ensure_size(ncol,idxmap%loc_to_glob,info,pad=-1_psb_lpk_,addsz=laddsz) + call psb_ensure_size(ncol,idxmap%loc_to_glob,info,& + & pad=-1_psb_lpk_,addsz=laddsz) if (info /= psb_success_) then info=1 call psb_errpush(psb_err_from_subroutine_ai_,name,& diff --git a/base/modules/psb_const_mod.F90 b/base/modules/psb_const_mod.F90 index c2277da0..85833a49 100644 --- a/base/modules/psb_const_mod.F90 +++ b/base/modules/psb_const_mod.F90 @@ -204,6 +204,9 @@ module psb_const_mod integer(psb_ipk_), parameter :: psb_spmat_null_=0, psb_spmat_bld_=1 integer(psb_ipk_), parameter :: psb_spmat_asb_=2, psb_spmat_upd_=4 + integer(psb_ipk_), parameter :: psb_matbld_noremote_=0, psb_matbld_remote_=1 + + integer(psb_ipk_), parameter :: psb_ireg_flgs_=10, psb_ip2_=0 integer(psb_ipk_), parameter :: psb_iflag_=2, psb_ichk_=3 integer(psb_ipk_), parameter :: psb_nnzt_=4, psb_zero_=5,psb_ipc_=6 diff --git a/base/modules/serial/psb_c_mat_mod.F90 b/base/modules/serial/psb_c_mat_mod.F90 index 75699249..2ba77fcf 100644 --- a/base/modules/serial/psb_c_mat_mod.F90 +++ b/base/modules/serial/psb_c_mat_mod.F90 @@ -85,6 +85,8 @@ module psb_c_mat_mod type :: psb_cspmat_type class(psb_c_base_sparse_mat), allocatable :: a + integer(psb_ipk_) :: remote_build=psb_matbld_noremote_ + class(psb_lc_base_sparse_mat), allocatable :: rmta contains ! Getters diff --git a/base/modules/serial/psb_d_mat_mod.F90 b/base/modules/serial/psb_d_mat_mod.F90 index ff51d1cb..d03818e7 100644 --- a/base/modules/serial/psb_d_mat_mod.F90 +++ b/base/modules/serial/psb_d_mat_mod.F90 @@ -85,6 +85,8 @@ module psb_d_mat_mod type :: psb_dspmat_type class(psb_d_base_sparse_mat), allocatable :: a + integer(psb_ipk_) :: remote_build=psb_matbld_noremote_ + class(psb_ld_base_sparse_mat), allocatable :: rmta contains ! Getters diff --git a/base/modules/serial/psb_s_mat_mod.F90 b/base/modules/serial/psb_s_mat_mod.F90 index 849c64c3..d748667a 100644 --- a/base/modules/serial/psb_s_mat_mod.F90 +++ b/base/modules/serial/psb_s_mat_mod.F90 @@ -85,6 +85,8 @@ module psb_s_mat_mod type :: psb_sspmat_type class(psb_s_base_sparse_mat), allocatable :: a + integer(psb_ipk_) :: remote_build=psb_matbld_noremote_ + class(psb_ls_base_sparse_mat), allocatable :: rmta contains ! Getters diff --git a/base/modules/serial/psb_z_mat_mod.F90 b/base/modules/serial/psb_z_mat_mod.F90 index fc16ca80..b81a5126 100644 --- a/base/modules/serial/psb_z_mat_mod.F90 +++ b/base/modules/serial/psb_z_mat_mod.F90 @@ -85,6 +85,8 @@ module psb_z_mat_mod type :: psb_zspmat_type class(psb_z_base_sparse_mat), allocatable :: a + integer(psb_ipk_) :: remote_build=psb_matbld_noremote_ + class(psb_lz_base_sparse_mat), allocatable :: rmta contains ! Getters diff --git a/base/modules/tools/psb_c_tools_mod.F90 b/base/modules/tools/psb_c_tools_mod.F90 index 378e146b..91a629e1 100644 --- a/base/modules/tools/psb_c_tools_mod.F90 +++ b/base/modules/tools/psb_c_tools_mod.F90 @@ -239,13 +239,13 @@ Module psb_c_tools_mod interface psb_spall - subroutine psb_cspalloc(a, desc_a, info, nnz) + subroutine psb_cspalloc(a, desc_a, info, nnz, bldmode) import implicit none type(psb_desc_type), intent(in) :: desc_a type(psb_cspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode end subroutine psb_cspalloc end interface diff --git a/base/modules/tools/psb_d_tools_mod.F90 b/base/modules/tools/psb_d_tools_mod.F90 index 81c75ece..f76a31eb 100644 --- a/base/modules/tools/psb_d_tools_mod.F90 +++ b/base/modules/tools/psb_d_tools_mod.F90 @@ -239,13 +239,13 @@ Module psb_d_tools_mod interface psb_spall - subroutine psb_dspalloc(a, desc_a, info, nnz) + subroutine psb_dspalloc(a, desc_a, info, nnz, bldmode) import implicit none type(psb_desc_type), intent(in) :: desc_a type(psb_dspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode end subroutine psb_dspalloc end interface diff --git a/base/modules/tools/psb_s_tools_mod.F90 b/base/modules/tools/psb_s_tools_mod.F90 index fa82a53e..aa2e6990 100644 --- a/base/modules/tools/psb_s_tools_mod.F90 +++ b/base/modules/tools/psb_s_tools_mod.F90 @@ -239,13 +239,13 @@ Module psb_s_tools_mod interface psb_spall - subroutine psb_sspalloc(a, desc_a, info, nnz) + subroutine psb_sspalloc(a, desc_a, info, nnz, bldmode) import implicit none type(psb_desc_type), intent(in) :: desc_a type(psb_sspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode end subroutine psb_sspalloc end interface diff --git a/base/modules/tools/psb_z_tools_mod.F90 b/base/modules/tools/psb_z_tools_mod.F90 index 233f2c20..bdb91e98 100644 --- a/base/modules/tools/psb_z_tools_mod.F90 +++ b/base/modules/tools/psb_z_tools_mod.F90 @@ -239,13 +239,13 @@ Module psb_z_tools_mod interface psb_spall - subroutine psb_zspalloc(a, desc_a, info, nnz) + subroutine psb_zspalloc(a, desc_a, info, nnz, bldmode) import implicit none type(psb_desc_type), intent(in) :: desc_a type(psb_zspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode end subroutine psb_zspalloc end interface diff --git a/base/serial/impl/psb_c_coo_impl.F90 b/base/serial/impl/psb_c_coo_impl.F90 index bd03b1b8..0127d00f 100644 --- a/base/serial/impl/psb_c_coo_impl.F90 +++ b/base/serial/impl/psb_c_coo_impl.F90 @@ -6121,7 +6121,6 @@ subroutine psb_lc_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info) if (nz == 0) return - nza = a%get_nzeros() isza = a%get_size() if (a%is_bld()) then diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90 index 088b012d..df5c4cd9 100644 --- a/base/serial/impl/psb_c_mat_impl.F90 +++ b/base/serial/impl/psb_c_mat_impl.F90 @@ -675,7 +675,12 @@ subroutine psb_c_free(a) call a%a%free() deallocate(a%a) endif - + if (allocated(a%rmta)) then + call a%rmta%free() + deallocate(a%rmta) + end if + a%remote_build = psb_matbld_noremote_ + end subroutine psb_c_free diff --git a/base/serial/impl/psb_d_coo_impl.F90 b/base/serial/impl/psb_d_coo_impl.F90 index d4d88027..7e835261 100644 --- a/base/serial/impl/psb_d_coo_impl.F90 +++ b/base/serial/impl/psb_d_coo_impl.F90 @@ -6121,7 +6121,6 @@ subroutine psb_ld_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info) if (nz == 0) return - nza = a%get_nzeros() isza = a%get_size() if (a%is_bld()) then diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90 index b21fa40f..2a6fb9a5 100644 --- a/base/serial/impl/psb_d_mat_impl.F90 +++ b/base/serial/impl/psb_d_mat_impl.F90 @@ -675,7 +675,12 @@ subroutine psb_d_free(a) call a%a%free() deallocate(a%a) endif - + if (allocated(a%rmta)) then + call a%rmta%free() + deallocate(a%rmta) + end if + a%remote_build = psb_matbld_noremote_ + end subroutine psb_d_free diff --git a/base/serial/impl/psb_s_coo_impl.F90 b/base/serial/impl/psb_s_coo_impl.F90 index 5ae9dd96..67725ffb 100644 --- a/base/serial/impl/psb_s_coo_impl.F90 +++ b/base/serial/impl/psb_s_coo_impl.F90 @@ -6121,7 +6121,6 @@ subroutine psb_ls_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info) if (nz == 0) return - nza = a%get_nzeros() isza = a%get_size() if (a%is_bld()) then diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90 index c624dc2f..ce7ce653 100644 --- a/base/serial/impl/psb_s_mat_impl.F90 +++ b/base/serial/impl/psb_s_mat_impl.F90 @@ -675,7 +675,12 @@ subroutine psb_s_free(a) call a%a%free() deallocate(a%a) endif - + if (allocated(a%rmta)) then + call a%rmta%free() + deallocate(a%rmta) + end if + a%remote_build = psb_matbld_noremote_ + end subroutine psb_s_free diff --git a/base/serial/impl/psb_z_coo_impl.F90 b/base/serial/impl/psb_z_coo_impl.F90 index 49e6c7fe..158f1ffe 100644 --- a/base/serial/impl/psb_z_coo_impl.F90 +++ b/base/serial/impl/psb_z_coo_impl.F90 @@ -6121,7 +6121,6 @@ subroutine psb_lz_coo_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info) if (nz == 0) return - nza = a%get_nzeros() isza = a%get_size() if (a%is_bld()) then diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90 index e008dfdb..2cebf9e7 100644 --- a/base/serial/impl/psb_z_mat_impl.F90 +++ b/base/serial/impl/psb_z_mat_impl.F90 @@ -675,7 +675,12 @@ subroutine psb_z_free(a) call a%a%free() deallocate(a%a) endif - + if (allocated(a%rmta)) then + call a%rmta%free() + deallocate(a%rmta) + end if + a%remote_build = psb_matbld_noremote_ + end subroutine psb_z_free diff --git a/base/tools/psb_cspalloc.f90 b/base/tools/psb_cspalloc.f90 index 69b9e1c2..d5a56609 100644 --- a/base/tools/psb_cspalloc.f90 +++ b/base/tools/psb_cspalloc.f90 @@ -41,7 +41,7 @@ ! nnz - integer(optional). The number of nonzeroes in the matrix. ! (local, user estimate) ! -subroutine psb_cspalloc(a, desc_a, info, nnz) +subroutine psb_cspalloc(a, desc_a, info, nnz, bldmode) use psb_base_mod, psb_protect_name => psb_cspalloc implicit none @@ -49,13 +49,13 @@ subroutine psb_cspalloc(a, desc_a, info, nnz) type(psb_desc_type), intent(in) :: desc_a type(psb_cspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode !locals type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, err_act - integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype - integer(psb_lpk_) :: m, n + integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype, bldmode_ + integer(psb_lpk_) :: m, n, nnzrmt_ integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -96,7 +96,12 @@ subroutine psb_cspalloc(a, desc_a, info, nnz) else nnz_ = max(1,5*loc_row) endif - + if (present(bldmode)) then + bldmode_ = bldmode + else + bldmode_ = psb_matbld_noremote_ + end if + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name), & & ':allocating size:',loc_row,loc_col,nnz_ @@ -109,6 +114,25 @@ subroutine psb_cspalloc(a, desc_a, info, nnz) goto 9999 end if + write(0,*) name,'Setting a%remote_build ',& + & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ + a%remote_build = bldmode_ + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + 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) + nnzrmt_ = max(100,(nnz_/100)) + call a%rmta%allocate(m,n,nnzrmt_) + + case default + write(0,*) name,'Invalid value for remote_build ' + a%remote_build = psb_matbld_noremote_ + end select + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),': ', & & desc_a%get_dectype(),psb_desc_bld_ diff --git a/base/tools/psb_cspasb.f90 b/base/tools/psb_cspasb.f90 index 96ed7fe7..0c2ade2d 100644 --- a/base/tools/psb_cspasb.f90 +++ b/base/tools/psb_cspasb.f90 @@ -111,7 +111,15 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) call a%set_ncols(n_col) end if - if (a%is_bld()) then + if (a%is_bld()) then + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + case (psb_matbld_remote_) + write(0,*) 'Need to implement data movement! ' + end select + call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then call a%asb(mold=mold) diff --git a/base/tools/psb_cspins.F90 b/base/tools/psb_cspins.F90 index 15ea556f..0084aa7d 100644 --- a/base/tools/psb_cspins.F90 +++ b/base/tools/psb_cspins.F90 @@ -70,6 +70,10 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) integer(psb_ipk_), parameter :: relocsz=200 logical :: rebuild_, local_ integer(psb_ipk_), allocatable :: ila(:),jla(:) + integer(psb_ipk_) :: i,k + integer(psb_lpk_) :: nnl + integer(psb_lpk_), allocatable :: lila(:),ljla(:) + complex(psb_spk_), allocatable :: lval(:) character(len=20) :: name info = psb_success_ @@ -147,6 +151,30 @@ 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) + 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) + + case default + write(0,*) name,' Ignoring wrong value for %remote_build' + end select + else info = psb_err_invalid_a_and_cd_state_ call psb_errpush(info,name) diff --git a/base/tools/psb_dspalloc.f90 b/base/tools/psb_dspalloc.f90 index 56ad6c93..088d5129 100644 --- a/base/tools/psb_dspalloc.f90 +++ b/base/tools/psb_dspalloc.f90 @@ -41,7 +41,7 @@ ! nnz - integer(optional). The number of nonzeroes in the matrix. ! (local, user estimate) ! -subroutine psb_dspalloc(a, desc_a, info, nnz) +subroutine psb_dspalloc(a, desc_a, info, nnz, bldmode) use psb_base_mod, psb_protect_name => psb_dspalloc implicit none @@ -49,13 +49,13 @@ subroutine psb_dspalloc(a, desc_a, info, nnz) type(psb_desc_type), intent(in) :: desc_a type(psb_dspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode !locals type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, err_act - integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype - integer(psb_lpk_) :: m, n + integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype, bldmode_ + integer(psb_lpk_) :: m, n, nnzrmt_ integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -96,7 +96,12 @@ subroutine psb_dspalloc(a, desc_a, info, nnz) else nnz_ = max(1,5*loc_row) endif - + if (present(bldmode)) then + bldmode_ = bldmode + else + bldmode_ = psb_matbld_noremote_ + end if + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name), & & ':allocating size:',loc_row,loc_col,nnz_ @@ -109,6 +114,25 @@ subroutine psb_dspalloc(a, desc_a, info, nnz) goto 9999 end if + write(0,*) name,'Setting a%remote_build ',& + & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ + a%remote_build = bldmode_ + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + 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) + nnzrmt_ = max(100,(nnz_/100)) + call a%rmta%allocate(m,n,nnzrmt_) + + case default + write(0,*) name,'Invalid value for remote_build ' + a%remote_build = psb_matbld_noremote_ + end select + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),': ', & & desc_a%get_dectype(),psb_desc_bld_ diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index 457553f7..293c5edd 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -111,7 +111,15 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) call a%set_ncols(n_col) end if - if (a%is_bld()) then + if (a%is_bld()) then + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + case (psb_matbld_remote_) + write(0,*) 'Need to implement data movement! ' + end select + call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then call a%asb(mold=mold) diff --git a/base/tools/psb_dspins.F90 b/base/tools/psb_dspins.F90 index 3e9ef0cc..1bd113dc 100644 --- a/base/tools/psb_dspins.F90 +++ b/base/tools/psb_dspins.F90 @@ -70,6 +70,10 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) integer(psb_ipk_), parameter :: relocsz=200 logical :: rebuild_, local_ integer(psb_ipk_), allocatable :: ila(:),jla(:) + integer(psb_ipk_) :: i,k + integer(psb_lpk_) :: nnl + integer(psb_lpk_), allocatable :: lila(:),ljla(:) + real(psb_dpk_), allocatable :: lval(:) character(len=20) :: name info = psb_success_ @@ -147,6 +151,30 @@ 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) + 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) + + case default + write(0,*) name,' Ignoring wrong value for %remote_build' + end select + else info = psb_err_invalid_a_and_cd_state_ call psb_errpush(info,name) diff --git a/base/tools/psb_sspalloc.f90 b/base/tools/psb_sspalloc.f90 index 15c3c538..1a418f2a 100644 --- a/base/tools/psb_sspalloc.f90 +++ b/base/tools/psb_sspalloc.f90 @@ -41,7 +41,7 @@ ! nnz - integer(optional). The number of nonzeroes in the matrix. ! (local, user estimate) ! -subroutine psb_sspalloc(a, desc_a, info, nnz) +subroutine psb_sspalloc(a, desc_a, info, nnz, bldmode) use psb_base_mod, psb_protect_name => psb_sspalloc implicit none @@ -49,13 +49,13 @@ subroutine psb_sspalloc(a, desc_a, info, nnz) type(psb_desc_type), intent(in) :: desc_a type(psb_sspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode !locals type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, err_act - integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype - integer(psb_lpk_) :: m, n + integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype, bldmode_ + integer(psb_lpk_) :: m, n, nnzrmt_ integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -96,7 +96,12 @@ subroutine psb_sspalloc(a, desc_a, info, nnz) else nnz_ = max(1,5*loc_row) endif - + if (present(bldmode)) then + bldmode_ = bldmode + else + bldmode_ = psb_matbld_noremote_ + end if + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name), & & ':allocating size:',loc_row,loc_col,nnz_ @@ -109,6 +114,25 @@ subroutine psb_sspalloc(a, desc_a, info, nnz) goto 9999 end if + write(0,*) name,'Setting a%remote_build ',& + & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ + a%remote_build = bldmode_ + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + 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) + nnzrmt_ = max(100,(nnz_/100)) + call a%rmta%allocate(m,n,nnzrmt_) + + case default + write(0,*) name,'Invalid value for remote_build ' + a%remote_build = psb_matbld_noremote_ + end select + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),': ', & & desc_a%get_dectype(),psb_desc_bld_ diff --git a/base/tools/psb_sspasb.f90 b/base/tools/psb_sspasb.f90 index f4e6169d..e27724b2 100644 --- a/base/tools/psb_sspasb.f90 +++ b/base/tools/psb_sspasb.f90 @@ -111,7 +111,15 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) call a%set_ncols(n_col) end if - if (a%is_bld()) then + if (a%is_bld()) then + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + case (psb_matbld_remote_) + write(0,*) 'Need to implement data movement! ' + end select + call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then call a%asb(mold=mold) diff --git a/base/tools/psb_sspins.F90 b/base/tools/psb_sspins.F90 index 1fd6eed0..2ae3d8df 100644 --- a/base/tools/psb_sspins.F90 +++ b/base/tools/psb_sspins.F90 @@ -70,6 +70,10 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) integer(psb_ipk_), parameter :: relocsz=200 logical :: rebuild_, local_ integer(psb_ipk_), allocatable :: ila(:),jla(:) + integer(psb_ipk_) :: i,k + integer(psb_lpk_) :: nnl + integer(psb_lpk_), allocatable :: lila(:),ljla(:) + real(psb_spk_), allocatable :: lval(:) character(len=20) :: name info = psb_success_ @@ -147,6 +151,30 @@ 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) + 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) + + case default + write(0,*) name,' Ignoring wrong value for %remote_build' + end select + else info = psb_err_invalid_a_and_cd_state_ call psb_errpush(info,name) diff --git a/base/tools/psb_zspalloc.f90 b/base/tools/psb_zspalloc.f90 index 16d48734..6d95f179 100644 --- a/base/tools/psb_zspalloc.f90 +++ b/base/tools/psb_zspalloc.f90 @@ -41,7 +41,7 @@ ! nnz - integer(optional). The number of nonzeroes in the matrix. ! (local, user estimate) ! -subroutine psb_zspalloc(a, desc_a, info, nnz) +subroutine psb_zspalloc(a, desc_a, info, nnz, bldmode) use psb_base_mod, psb_protect_name => psb_zspalloc implicit none @@ -49,13 +49,13 @@ subroutine psb_zspalloc(a, desc_a, info, nnz) type(psb_desc_type), intent(in) :: desc_a type(psb_zspmat_type), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: nnz + integer(psb_ipk_), optional, intent(in) :: nnz, bldmode !locals type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, err_act - integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype - integer(psb_lpk_) :: m, n + integer(psb_ipk_) :: loc_row,loc_col, nnz_, dectype, bldmode_ + integer(psb_lpk_) :: m, n, nnzrmt_ integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -96,7 +96,12 @@ subroutine psb_zspalloc(a, desc_a, info, nnz) else nnz_ = max(1,5*loc_row) endif - + if (present(bldmode)) then + bldmode_ = bldmode + else + bldmode_ = psb_matbld_noremote_ + end if + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name), & & ':allocating size:',loc_row,loc_col,nnz_ @@ -109,6 +114,25 @@ subroutine psb_zspalloc(a, desc_a, info, nnz) goto 9999 end if + write(0,*) name,'Setting a%remote_build ',& + & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ + a%remote_build = bldmode_ + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + 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) + nnzrmt_ = max(100,(nnz_/100)) + call a%rmta%allocate(m,n,nnzrmt_) + + case default + write(0,*) name,'Invalid value for remote_build ' + a%remote_build = psb_matbld_noremote_ + end select + if (debug_level >= psb_debug_ext_) & & write(debug_unit,*) me,' ',trim(name),': ', & & desc_a%get_dectype(),psb_desc_bld_ diff --git a/base/tools/psb_zspasb.f90 b/base/tools/psb_zspasb.f90 index b5966110..7021bbaa 100644 --- a/base/tools/psb_zspasb.f90 +++ b/base/tools/psb_zspasb.f90 @@ -111,7 +111,15 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) call a%set_ncols(n_col) end if - if (a%is_bld()) then + if (a%is_bld()) then + + select case(a%remote_build) + case (psb_matbld_noremote_) + ! nothing needed + case (psb_matbld_remote_) + write(0,*) 'Need to implement data movement! ' + end select + call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then call a%asb(mold=mold) diff --git a/base/tools/psb_zspins.F90 b/base/tools/psb_zspins.F90 index 525ed415..d76450e3 100644 --- a/base/tools/psb_zspins.F90 +++ b/base/tools/psb_zspins.F90 @@ -70,6 +70,10 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) integer(psb_ipk_), parameter :: relocsz=200 logical :: rebuild_, local_ integer(psb_ipk_), allocatable :: ila(:),jla(:) + integer(psb_ipk_) :: i,k + integer(psb_lpk_) :: nnl + integer(psb_lpk_), allocatable :: lila(:),ljla(:) + complex(psb_dpk_), allocatable :: lval(:) character(len=20) :: name info = psb_success_ @@ -147,6 +151,30 @@ 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) + 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) + + case default + write(0,*) name,' Ignoring wrong value for %remote_build' + end select + else info = psb_err_invalid_a_and_cd_state_ call psb_errpush(info,name) diff --git a/cbind/base/psb_c_serial_cbind_mod.F90 b/cbind/base/psb_c_serial_cbind_mod.F90 index dba41b67..b298d84a 100644 --- a/cbind/base/psb_c_serial_cbind_mod.F90 +++ b/cbind/base/psb_c_serial_cbind_mod.F90 @@ -204,38 +204,5 @@ contains end function psb_c_cvect_set_vect - function psb_c_g2l(cdh,gindex,cowned) bind(c) result(lindex) - use psb_base_mod - implicit none - - integer(psb_c_lpk_), value :: gindex - logical(c_bool), value :: cowned - type(psb_c_descriptor) :: cdh - integer(psb_c_ipk_) :: lindex - - type(psb_desc_type), pointer :: descp - integer(psb_ipk_) :: info, localindex, ixb, iam, np - logical :: owned - - ixb = psb_c_get_index_base() - owned = cowned - lindex = -1 - if (c_associated(cdh%item)) then - call c_f_pointer(cdh%item,descp) - else - return - end if - - call psb_info(descp%get_context(),iam,np) - if (ixb == 1) then - call descp%indxmap%g2l(gindex,localindex,info,owned=owned) - lindex = localindex - else - call descp%indxmap%g2l(gindex+(1-ixb),localindex,info,owned=owned) - lindex = localindex-(1-ixb) - endif - - end function psb_c_g2l - end module psb_c_serial_cbind_mod diff --git a/test/pargen/psb_d_pde3d.F90 b/test/pargen/psb_d_pde3d.F90 index def2ea37..2e3976ce 100644 --- a/test/pargen/psb_d_pde3d.F90 +++ b/test/pargen/psb_d_pde3d.F90 @@ -429,7 +429,8 @@ contains end select - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz, & + & bldmode=psb_matbld_remote_) ! define rhs from boundary conditions; also build initial guess if (info == psb_success_) call psb_geall(xv,desc_a,info) if (info == psb_success_) call psb_geall(bv,desc_a,info) diff --git a/test/pargen/psb_s_pde3d.F90 b/test/pargen/psb_s_pde3d.F90 index 619a7f75..a4788338 100644 --- a/test/pargen/psb_s_pde3d.F90 +++ b/test/pargen/psb_s_pde3d.F90 @@ -429,7 +429,8 @@ contains end select - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz, & + & bldmode=psb_matbld_remote_) ! define rhs from boundary conditions; also build initial guess if (info == psb_success_) call psb_geall(xv,desc_a,info) if (info == psb_success_) call psb_geall(bv,desc_a,info) From 0e676d2903c8f57ef6247a88ddf8db1826a1542e Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 1 Feb 2022 13:15:48 +0100 Subject: [PATCH 02/28] Infrastructure for remote builds --- base/tools/psb_cspasb.f90 | 2 +- base/tools/psb_cspins.F90 | 30 +++++++++++++------------- base/tools/psb_dspasb.f90 | 2 +- base/tools/psb_dspins.F90 | 31 ++++++++++++++------------- base/tools/psb_sspasb.f90 | 2 +- base/tools/psb_sspins.F90 | 30 +++++++++++++------------- base/tools/psb_zspasb.f90 | 2 +- base/tools/psb_zspins.F90 | 30 +++++++++++++------------- test/pargen/psb_d_pde3d.F90 | 42 ++++++++++++++++++++++++++++++++----- 9 files changed, 106 insertions(+), 65 deletions(-) diff --git a/base/tools/psb_cspasb.f90 b/base/tools/psb_cspasb.f90 index 0c2ade2d..747d1427 100644 --- a/base/tools/psb_cspasb.f90 +++ b/base/tools/psb_cspasb.f90 @@ -117,7 +117,7 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) 'Need to implement data movement! ' + write(0,*) me,' Size of rmta:',a%rmta%get_nzeros() end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_cspins.F90 b/base/tools/psb_cspins.F90 index 0084aa7d..2961c646 100644 --- a/base/tools/psb_cspins.F90 +++ b/base/tools/psb_cspins.F90 @@ -157,20 +157,22 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) ! Do nothing case (psb_matbld_remote_) nnl = count(ila(1:nz)<0) - 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) - + 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 diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index 293c5edd..14204ad9 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -117,7 +117,7 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) 'Need to implement data movement! ' + write(0,*) me,' Size of rmta:',a%rmta%get_nzeros() end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_dspins.F90 b/base/tools/psb_dspins.F90 index 1bd113dc..ac074487 100644 --- a/base/tools/psb_dspins.F90 +++ b/base/tools/psb_dspins.F90 @@ -156,21 +156,24 @@ 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) - 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) + 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 diff --git a/base/tools/psb_sspasb.f90 b/base/tools/psb_sspasb.f90 index e27724b2..4c019622 100644 --- a/base/tools/psb_sspasb.f90 +++ b/base/tools/psb_sspasb.f90 @@ -117,7 +117,7 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) 'Need to implement data movement! ' + write(0,*) me,' Size of rmta:',a%rmta%get_nzeros() end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_sspins.F90 b/base/tools/psb_sspins.F90 index 2ae3d8df..f8971ce2 100644 --- a/base/tools/psb_sspins.F90 +++ b/base/tools/psb_sspins.F90 @@ -157,20 +157,22 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) ! Do nothing case (psb_matbld_remote_) nnl = count(ila(1:nz)<0) - 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) - + 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 diff --git a/base/tools/psb_zspasb.f90 b/base/tools/psb_zspasb.f90 index 7021bbaa..6c9fc76f 100644 --- a/base/tools/psb_zspasb.f90 +++ b/base/tools/psb_zspasb.f90 @@ -117,7 +117,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) 'Need to implement data movement! ' + write(0,*) me,' Size of rmta:',a%rmta%get_nzeros() end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_zspins.F90 b/base/tools/psb_zspins.F90 index d76450e3..bfe61a17 100644 --- a/base/tools/psb_zspins.F90 +++ b/base/tools/psb_zspins.F90 @@ -157,20 +157,22 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) ! Do nothing case (psb_matbld_remote_) nnl = count(ila(1:nz)<0) - 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) - + 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 diff --git a/test/pargen/psb_d_pde3d.F90 b/test/pargen/psb_d_pde3d.F90 index 2e3976ce..bbef8645 100644 --- a/test/pargen/psb_d_pde3d.F90 +++ b/test/pargen/psb_d_pde3d.F90 @@ -208,7 +208,7 @@ contains type(psb_d_coo_sparse_mat) :: acoo type(psb_d_csr_sparse_mat) :: acsr real(psb_dpk_) :: 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,30 @@ 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 ! @@ -463,8 +487,9 @@ contains call psb_barrier(ctxt) t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) + do ii=1, mysz, nb + !ib = min(nb,nlr-ii+1) + ib = min(nb,mysz-ii+1) icoeff = 1 do k=1,ib i=ii+k-1 @@ -616,8 +641,15 @@ contains return end subroutine psb_d_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+1)) & + & .or.(j=bndy(iamy+1)) & + & .or.(k=bndz(iamz+1)) + end function outside end module psb_d_pde3d_mod program psb_d_pde3d From 6d0b26ecf13e8c468bb57526d2cb739c0511aeeb Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 3 Feb 2022 13:25:51 +0100 Subject: [PATCH 03/28] Further changes for remote build, new method --- base/internals/psi_graph_fnd_owner.F90 | 12 +- base/internals/psi_indx_map_fnd_owner.F90 | 10 +- base/modules/desc/psb_indx_map_mod.f90 | 3 +- base/modules/serial/psb_c_mat_mod.F90 | 2 +- base/modules/serial/psb_d_mat_mod.F90 | 2 +- base/modules/serial/psb_s_mat_mod.F90 | 2 +- base/modules/serial/psb_z_mat_mod.F90 | 2 +- base/modules/tools/psb_c_tools_mod.F90 | 13 +- base/modules/tools/psb_d_tools_mod.F90 | 13 +- base/modules/tools/psb_s_tools_mod.F90 | 13 +- base/modules/tools/psb_z_tools_mod.F90 | 13 +- base/tools/Makefile | 3 +- base/tools/psb_c_remote_mat.F90 | 259 ++++++++++++++++++++++ base/tools/psb_cspalloc.f90 | 2 +- base/tools/psb_cspasb.f90 | 12 +- base/tools/psb_cspins.F90 | 32 ++- base/tools/psb_d_remote_mat.F90 | 259 ++++++++++++++++++++++ base/tools/psb_dspalloc.f90 | 2 +- base/tools/psb_dspasb.f90 | 12 +- base/tools/psb_dspins.F90 | 33 ++- base/tools/psb_s_remote_mat.F90 | 259 ++++++++++++++++++++++ base/tools/psb_sspalloc.f90 | 2 +- base/tools/psb_sspasb.f90 | 12 +- base/tools/psb_sspins.F90 | 32 ++- base/tools/psb_z_remote_mat.F90 | 259 ++++++++++++++++++++++ base/tools/psb_zspalloc.f90 | 2 +- base/tools/psb_zspasb.f90 | 12 +- base/tools/psb_zspins.F90 | 32 ++- test/pargen/psb_d_pde3d.F90 | 141 +++++++++--- test/pargen/psb_s_pde3d.F90 | 127 ++++++++++- 30 files changed, 1499 insertions(+), 78 deletions(-) create mode 100644 base/tools/psb_c_remote_mat.F90 create mode 100644 base/tools/psb_d_remote_mat.F90 create mode 100644 base/tools/psb_s_remote_mat.F90 create mode 100644 base/tools/psb_z_remote_mat.F90 diff --git a/base/internals/psi_graph_fnd_owner.F90 b/base/internals/psi_graph_fnd_owner.F90 index 966af403..e3e7ec4f 100644 --- a/base/internals/psi_graph_fnd_owner.F90 +++ b/base/internals/psi_graph_fnd_owner.F90 @@ -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) ! diff --git a/base/internals/psi_indx_map_fnd_owner.F90 b/base/internals/psi_indx_map_fnd_owner.F90 index 0ecade3d..3efcbe59 100644 --- a/base/internals/psi_indx_map_fnd_owner.F90 +++ b/base/internals/psi_indx_map_fnd_owner.F90 @@ -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 diff --git a/base/modules/desc/psb_indx_map_mod.f90 b/base/modules/desc/psb_indx_map_mod.f90 index 79526d59..d4ce78f1 100644 --- a/base/modules/desc/psb_indx_map_mod.f90 +++ b/base/modules/desc/psb_indx_map_mod.f90 @@ -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 diff --git a/base/modules/serial/psb_c_mat_mod.F90 b/base/modules/serial/psb_c_mat_mod.F90 index 2ba77fcf..e3e8ab69 100644 --- a/base/modules/serial/psb_c_mat_mod.F90 +++ b/base/modules/serial/psb_c_mat_mod.F90 @@ -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 diff --git a/base/modules/serial/psb_d_mat_mod.F90 b/base/modules/serial/psb_d_mat_mod.F90 index d03818e7..a072197e 100644 --- a/base/modules/serial/psb_d_mat_mod.F90 +++ b/base/modules/serial/psb_d_mat_mod.F90 @@ -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 diff --git a/base/modules/serial/psb_s_mat_mod.F90 b/base/modules/serial/psb_s_mat_mod.F90 index d748667a..20bf3249 100644 --- a/base/modules/serial/psb_s_mat_mod.F90 +++ b/base/modules/serial/psb_s_mat_mod.F90 @@ -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 diff --git a/base/modules/serial/psb_z_mat_mod.F90 b/base/modules/serial/psb_z_mat_mod.F90 index b81a5126..1ca6941b 100644 --- a/base/modules/serial/psb_z_mat_mod.F90 +++ b/base/modules/serial/psb_z_mat_mod.F90 @@ -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 diff --git a/base/modules/tools/psb_c_tools_mod.F90 b/base/modules/tools/psb_c_tools_mod.F90 index 91a629e1..1aa2bcbb 100644 --- a/base/modules/tools/psb_c_tools_mod.F90 +++ b/base/modules/tools/psb_c_tools_mod.F90 @@ -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 diff --git a/base/modules/tools/psb_d_tools_mod.F90 b/base/modules/tools/psb_d_tools_mod.F90 index f76a31eb..b73fcb74 100644 --- a/base/modules/tools/psb_d_tools_mod.F90 +++ b/base/modules/tools/psb_d_tools_mod.F90 @@ -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 diff --git a/base/modules/tools/psb_s_tools_mod.F90 b/base/modules/tools/psb_s_tools_mod.F90 index aa2e6990..dfa18e92 100644 --- a/base/modules/tools/psb_s_tools_mod.F90 +++ b/base/modules/tools/psb_s_tools_mod.F90 @@ -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 diff --git a/base/modules/tools/psb_z_tools_mod.F90 b/base/modules/tools/psb_z_tools_mod.F90 index bdb91e98..b76893c7 100644 --- a/base/modules/tools/psb_z_tools_mod.F90 +++ b/base/modules/tools/psb_z_tools_mod.F90 @@ -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 diff --git a/base/tools/Makefile b/base/tools/Makefile index 2c160616..ce13caed 100644 --- a/base/tools/Makefile +++ b/base/tools/Makefile @@ -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=.. diff --git a/base/tools/psb_c_remote_mat.F90 b/base/tools/psb_c_remote_mat.F90 new file mode 100644 index 00000000..acd2ba6e --- /dev/null +++ b/base/tools/psb_c_remote_mat.F90 @@ -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 + diff --git a/base/tools/psb_cspalloc.f90 b/base/tools/psb_cspalloc.f90 index d5a56609..379cb935 100644 --- a/base/tools/psb_cspalloc.f90 +++ b/base/tools/psb_cspalloc.f90 @@ -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_) diff --git a/base/tools/psb_cspasb.f90 b/base/tools/psb_cspasb.f90 index 747d1427..4904d336 100644 --- a/base/tools/psb_cspasb.f90 +++ b/base/tools/psb_cspasb.f90 @@ -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) diff --git a/base/tools/psb_cspins.F90 b/base/tools/psb_cspins.F90 index 2961c646..6c050b0c 100644 --- a/base/tools/psb_cspins.F90 +++ b/base/tools/psb_cspins.F90 @@ -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) diff --git a/base/tools/psb_d_remote_mat.F90 b/base/tools/psb_d_remote_mat.F90 new file mode 100644 index 00000000..3e35970a --- /dev/null +++ b/base/tools/psb_d_remote_mat.F90 @@ -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 + diff --git a/base/tools/psb_dspalloc.f90 b/base/tools/psb_dspalloc.f90 index 088d5129..3d945274 100644 --- a/base/tools/psb_dspalloc.f90 +++ b/base/tools/psb_dspalloc.f90 @@ -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_) diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index 14204ad9..bc66eb11 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -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) diff --git a/base/tools/psb_dspins.F90 b/base/tools/psb_dspins.F90 index ac074487..06276fa8 100644 --- a/base/tools/psb_dspins.F90 +++ b/base/tools/psb_dspins.F90 @@ -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) diff --git a/base/tools/psb_s_remote_mat.F90 b/base/tools/psb_s_remote_mat.F90 new file mode 100644 index 00000000..0a78d2b6 --- /dev/null +++ b/base/tools/psb_s_remote_mat.F90 @@ -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 + diff --git a/base/tools/psb_sspalloc.f90 b/base/tools/psb_sspalloc.f90 index 1a418f2a..2d79b0ba 100644 --- a/base/tools/psb_sspalloc.f90 +++ b/base/tools/psb_sspalloc.f90 @@ -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_) diff --git a/base/tools/psb_sspasb.f90 b/base/tools/psb_sspasb.f90 index 4c019622..fa9a2a49 100644 --- a/base/tools/psb_sspasb.f90 +++ b/base/tools/psb_sspasb.f90 @@ -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) diff --git a/base/tools/psb_sspins.F90 b/base/tools/psb_sspins.F90 index f8971ce2..c4abb9f6 100644 --- a/base/tools/psb_sspins.F90 +++ b/base/tools/psb_sspins.F90 @@ -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) diff --git a/base/tools/psb_z_remote_mat.F90 b/base/tools/psb_z_remote_mat.F90 new file mode 100644 index 00000000..1acc57b2 --- /dev/null +++ b/base/tools/psb_z_remote_mat.F90 @@ -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 + diff --git a/base/tools/psb_zspalloc.f90 b/base/tools/psb_zspalloc.f90 index 6d95f179..1d76d66c 100644 --- a/base/tools/psb_zspalloc.f90 +++ b/base/tools/psb_zspalloc.f90 @@ -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_) diff --git a/base/tools/psb_zspasb.f90 b/base/tools/psb_zspasb.f90 index 6c9fc76f..67353739 100644 --- a/base/tools/psb_zspasb.f90 +++ b/base/tools/psb_zspasb.f90 @@ -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) diff --git a/base/tools/psb_zspins.F90 b/base/tools/psb_zspins.F90 index bfe61a17..2a654e14 100644 --- a/base/tools/psb_zspins.F90 +++ b/base/tools/psb_zspins.F90 @@ -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) diff --git a/test/pargen/psb_d_pde3d.F90 b/test/pargen/psb_d_pde3d.F90 index bbef8645..483e6d9c 100644 --- a/test/pargen/psb_d_pde3d.F90 +++ b/test/pargen/psb_d_pde3d.F90 @@ -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 diff --git a/test/pargen/psb_s_pde3d.F90 b/test/pargen/psb_s_pde3d.F90 index a4788338..b27f8818 100644 --- a/test/pargen/psb_s_pde3d.F90 +++ b/test/pargen/psb_s_pde3d.F90 @@ -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+1)) & + & .or.(j=bndy(iamy+1)) & + & .or.(k=bndz(iamz+1)) + end function outside end module psb_s_pde3d_mod program psb_s_pde3d From 0b19adab3c242d3d85614e0d2a27fafa8f283a17 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 4 Feb 2022 19:12:32 +0100 Subject: [PATCH 04/28] Round of fixes for matrix remote build. Something still wrong. --- base/tools/psb_c_remote_mat.F90 | 147 ++++++++++++++++++-------------- base/tools/psb_cspalloc.f90 | 8 +- base/tools/psb_cspasb.f90 | 36 +++++--- base/tools/psb_cspins.F90 | 6 +- base/tools/psb_d_remote_mat.F90 | 147 ++++++++++++++++++-------------- base/tools/psb_dspalloc.f90 | 8 +- base/tools/psb_dspasb.f90 | 36 +++++--- base/tools/psb_dspins.F90 | 6 +- base/tools/psb_s_remote_mat.F90 | 147 ++++++++++++++++++-------------- base/tools/psb_sspalloc.f90 | 8 +- base/tools/psb_sspasb.f90 | 36 +++++--- base/tools/psb_sspins.F90 | 6 +- base/tools/psb_z_remote_mat.F90 | 147 ++++++++++++++++++-------------- base/tools/psb_zspalloc.f90 | 8 +- base/tools/psb_zspasb.f90 | 36 +++++--- base/tools/psb_zspins.F90 | 6 +- test/pargen/psb_d_pde3d.F90 | 4 +- test/pargen/psb_s_pde3d.F90 | 4 +- 18 files changed, 464 insertions(+), 332 deletions(-) diff --git a/base/tools/psb_c_remote_mat.F90 b/base/tools/psb_c_remote_mat.F90 index acd2ba6e..fcccd26d 100644 --- a/base/tools/psb_c_remote_mat.F90 +++ b/base/tools/psb_c_remote_mat.F90 @@ -98,20 +98,19 @@ Subroutine psb_lc_remote_mat(a,desc_a,b,info) integer(psb_ipk_) :: nnp, nrcvs, nsnds integer(psb_mpk_) :: icomm, minfo integer(psb_mpk_), allocatable :: brvindx(:), & - & rvsz(:), bsdindx(:),sdsz(:) + & rvsz(:), bsdindx(:),sdsz(:), sdsi(:), rvsi(:) 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(:) + integer(psb_ipk_), allocatable :: 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' + name='psb_c_remote_mat' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then info = psb_err_internal_error_ ; goto 9999 @@ -129,9 +128,9 @@ Subroutine psb_lc_remote_mat(a,desc_a,b,info) call b%free() - - Allocate(brvindx(np+1),& - & rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info) + + Allocate(rvsz(np),sdsz(np),sdsi(np),rvsi(np),brvindx(np+1),& + & bsdindx(np+1), acoo,stat=info) if (info /= psb_success_) then info=psb_err_alloc_dealloc_ @@ -140,71 +139,94 @@ Subroutine psb_lc_remote_mat(a,desc_a,b,info) 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) + 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' - 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 - + 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) - - - lnnz = max(iszr,iszs,lone) - lnc = a%get_ncols() - call acoo%allocate(lnr,lnc,lnnz) + icomm = desc_a%get_mpic() + sdsz(:)=0 + rvsz(:)=0 + sdsi(:)=0 + rvsi(:)=0 + ipx = 1 + brvindx(:) = 0 + bsdindx(:) = 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 + Do proc=0,np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + iszs = sum(sdsz) + iszr = sum(rvsz) + call acoo%allocate(desc_a%get_global_rows(),desc_a%get_global_cols(),iszr) + if (psb_errstatus_fatal()) then + write(0,*) 'Error from acoo%allocate ' + info = 4010 + goto 9999 + end if if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& & ' Send:',sdsz(:),' Receive:',rvsz(:) - + !write(debug_unit,*) me,' ',trim(name),': ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info) + !write(debug_unit,*) me,' ',trim(name),' iasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + !write(debug_unit,*) me,' ',trim(name),' jasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + !write(debug_unit,*) me,' ',trim(name),' 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 + do k=1, nz + proc = iprc(k) + sdsi(proc+1) = sdsi(proc+1) + 1 + !rvsi(proc) = rvsi(proc) + 1 + iasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ia(k) + jasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ja(k) + valsnd(bsdindx(proc+1)+sdsi(proc+1)) = a%val(k) + end do + do proc=0,np-1 + if (sdsi(proc+1) /= sdsz(proc+1)) & + & write(0,*) me,name,'Send mismacth ',sdsi(proc+1),sdsz(proc+1) + end do select case(psb_get_sp_a2av_alg()) case(psb_sp_a2av_smpl_triad_) @@ -218,7 +240,7 @@ Subroutine psb_lc_remote_mat(a,desc_a,b,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) & @@ -239,10 +261,9 @@ Subroutine psb_lc_remote_mat(a,desc_a,b,info) call psb_errpush(info,name,a_err='alltoallv') goto 9999 end if - + call acoo%set_nzeros(iszr) call acoo%mv_to_coo(b,info) - - + Deallocate(brvindx,bsdindx,rvsz,sdsz,& & iasnd,jasnd,valsnd,stat=info) if (debug_level >= psb_debug_outer_)& diff --git a/base/tools/psb_cspalloc.f90 b/base/tools/psb_cspalloc.f90 index 379cb935..286bb4a0 100644 --- a/base/tools/psb_cspalloc.f90 +++ b/base/tools/psb_cspalloc.f90 @@ -114,16 +114,16 @@ subroutine psb_cspalloc(a, desc_a, info, nnz, bldmode) goto 9999 end if - write(0,*) name,'Setting a%remote_build ',& - & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ +!!$ write(0,*) name,'Setting a%remote_build ',& +!!$ & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ a%remote_build = bldmode_ select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed - write(0,*) name,' matbld_noremote_ nothing needed' + !write(0,*) name,' matbld_noremote_ nothing needed' case (psb_matbld_remote_) - write(0,*) name,' matbld_remote_ start ' + !write(0,*) name,' matbld_remote_ start ' allocate(a%rmta) nnzrmt_ = max(100,(nnz_/100)) call a%rmta%allocate(m,n,nnzrmt_) diff --git a/base/tools/psb_cspasb.f90 b/base/tools/psb_cspasb.f90 index 4904d336..99c8990e 100644 --- a/base/tools/psb_cspasb.f90 +++ b/base/tools/psb_cspasb.f90 @@ -106,26 +106,38 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) ! First case: we come from a fresh build. ! - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - call a%set_nrows(n_row) - call a%set_ncols(n_col) - end if - - if (a%is_bld()) then - select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) me,name,' 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 - + integer(psb_ipk_), allocatable :: ila(:), jla(:) + integer(psb_ipk_) :: nz, nzt,k call psb_remote_mat(a%rmta,desc_a,a_add,info) - - end block + nz = a_add%get_nzeros() +!!$ write(0,*) me,name,' Nz to be added',nz + nzt = nz + call psb_sum(ctxt,nzt) + if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nz > 0) then + ! + ! Should we check for new indices here? + ! + call psb_realloc(nz,ila,info) + call psb_realloc(nz,jla,info) + call desc_a%indxmap%g2l(a_add%ia(1:nz),ila(1:nz),info,owned=.true.) + if (info == 0) call desc_a%indxmap%g2l_ins(a_add%ja(1:nz),jla(1:nz),info) + !write(0,*) me,name,' Check before insert',a%get_nzeros() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) + !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz + end if + if (nzt > 0) call psb_cdasb(desc_a,info) + end block end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_cspins.F90 b/base/tools/psb_cspins.F90 index 6c050b0c..553e3e3d 100644 --- a/base/tools/psb_cspins.F90 +++ b/base/tools/psb_cspins.F90 @@ -164,9 +164,9 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) do i=1,nz if (ila(i)<0) then k=k+1 - lila(k) = ia(k) - ljla(k) = ja(k) - lval(k) = val(k) + lila(k) = ia(i) + ljla(k) = ja(i) + lval(k) = val(i) end if end do if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl diff --git a/base/tools/psb_d_remote_mat.F90 b/base/tools/psb_d_remote_mat.F90 index 3e35970a..290a6f3d 100644 --- a/base/tools/psb_d_remote_mat.F90 +++ b/base/tools/psb_d_remote_mat.F90 @@ -98,20 +98,19 @@ Subroutine psb_ld_remote_mat(a,desc_a,b,info) integer(psb_ipk_) :: nnp, nrcvs, nsnds integer(psb_mpk_) :: icomm, minfo integer(psb_mpk_), allocatable :: brvindx(:), & - & rvsz(:), bsdindx(:),sdsz(:) + & rvsz(:), bsdindx(:),sdsz(:), sdsi(:), rvsi(:) 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(:) + integer(psb_ipk_), allocatable :: 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' + name='psb_d_remote_mat' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then info = psb_err_internal_error_ ; goto 9999 @@ -129,9 +128,9 @@ Subroutine psb_ld_remote_mat(a,desc_a,b,info) call b%free() - - Allocate(brvindx(np+1),& - & rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info) + + Allocate(rvsz(np),sdsz(np),sdsi(np),rvsi(np),brvindx(np+1),& + & bsdindx(np+1), acoo,stat=info) if (info /= psb_success_) then info=psb_err_alloc_dealloc_ @@ -140,71 +139,94 @@ Subroutine psb_ld_remote_mat(a,desc_a,b,info) 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) + 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' - 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 - + 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) - - - lnnz = max(iszr,iszs,lone) - lnc = a%get_ncols() - call acoo%allocate(lnr,lnc,lnnz) + icomm = desc_a%get_mpic() + sdsz(:)=0 + rvsz(:)=0 + sdsi(:)=0 + rvsi(:)=0 + ipx = 1 + brvindx(:) = 0 + bsdindx(:) = 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 + Do proc=0,np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + iszs = sum(sdsz) + iszr = sum(rvsz) + call acoo%allocate(desc_a%get_global_rows(),desc_a%get_global_cols(),iszr) + if (psb_errstatus_fatal()) then + write(0,*) 'Error from acoo%allocate ' + info = 4010 + goto 9999 + end if if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& & ' Send:',sdsz(:),' Receive:',rvsz(:) - + !write(debug_unit,*) me,' ',trim(name),': ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info) + !write(debug_unit,*) me,' ',trim(name),' iasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + !write(debug_unit,*) me,' ',trim(name),' jasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + !write(debug_unit,*) me,' ',trim(name),' 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 + do k=1, nz + proc = iprc(k) + sdsi(proc+1) = sdsi(proc+1) + 1 + !rvsi(proc) = rvsi(proc) + 1 + iasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ia(k) + jasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ja(k) + valsnd(bsdindx(proc+1)+sdsi(proc+1)) = a%val(k) + end do + do proc=0,np-1 + if (sdsi(proc+1) /= sdsz(proc+1)) & + & write(0,*) me,name,'Send mismacth ',sdsi(proc+1),sdsz(proc+1) + end do select case(psb_get_sp_a2av_alg()) case(psb_sp_a2av_smpl_triad_) @@ -218,7 +240,7 @@ Subroutine psb_ld_remote_mat(a,desc_a,b,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) & @@ -239,10 +261,9 @@ Subroutine psb_ld_remote_mat(a,desc_a,b,info) call psb_errpush(info,name,a_err='alltoallv') goto 9999 end if - + call acoo%set_nzeros(iszr) call acoo%mv_to_coo(b,info) - - + Deallocate(brvindx,bsdindx,rvsz,sdsz,& & iasnd,jasnd,valsnd,stat=info) if (debug_level >= psb_debug_outer_)& diff --git a/base/tools/psb_dspalloc.f90 b/base/tools/psb_dspalloc.f90 index 3d945274..781a3abf 100644 --- a/base/tools/psb_dspalloc.f90 +++ b/base/tools/psb_dspalloc.f90 @@ -114,16 +114,16 @@ subroutine psb_dspalloc(a, desc_a, info, nnz, bldmode) goto 9999 end if - write(0,*) name,'Setting a%remote_build ',& - & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ +!!$ write(0,*) name,'Setting a%remote_build ',& +!!$ & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ a%remote_build = bldmode_ select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed - write(0,*) name,' matbld_noremote_ nothing needed' + !write(0,*) name,' matbld_noremote_ nothing needed' case (psb_matbld_remote_) - write(0,*) name,' matbld_remote_ start ' + !write(0,*) name,' matbld_remote_ start ' allocate(a%rmta) nnzrmt_ = max(100,(nnz_/100)) call a%rmta%allocate(m,n,nnzrmt_) diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index bc66eb11..c9d3cd08 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -106,26 +106,38 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) ! First case: we come from a fresh build. ! - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - call a%set_nrows(n_row) - call a%set_ncols(n_col) - end if - - if (a%is_bld()) then - select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) me,name,' 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 - + integer(psb_ipk_), allocatable :: ila(:), jla(:) + integer(psb_ipk_) :: nz, nzt,k call psb_remote_mat(a%rmta,desc_a,a_add,info) - - end block + nz = a_add%get_nzeros() +!!$ write(0,*) me,name,' Nz to be added',nz + nzt = nz + call psb_sum(ctxt,nzt) + if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nz > 0) then + ! + ! Should we check for new indices here? + ! + call psb_realloc(nz,ila,info) + call psb_realloc(nz,jla,info) + call desc_a%indxmap%g2l(a_add%ia(1:nz),ila(1:nz),info,owned=.true.) + if (info == 0) call desc_a%indxmap%g2l_ins(a_add%ja(1:nz),jla(1:nz),info) + !write(0,*) me,name,' Check before insert',a%get_nzeros() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) + !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz + end if + if (nzt > 0) call psb_cdasb(desc_a,info) + end block end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_dspins.F90 b/base/tools/psb_dspins.F90 index 06276fa8..064b816c 100644 --- a/base/tools/psb_dspins.F90 +++ b/base/tools/psb_dspins.F90 @@ -164,9 +164,9 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) do i=1,nz if (ila(i)<0) then k=k+1 - lila(k) = ia(k) - ljla(k) = ja(k) - lval(k) = val(k) + lila(k) = ia(i) + ljla(k) = ja(i) + lval(k) = val(i) end if end do if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl diff --git a/base/tools/psb_s_remote_mat.F90 b/base/tools/psb_s_remote_mat.F90 index 0a78d2b6..26b12652 100644 --- a/base/tools/psb_s_remote_mat.F90 +++ b/base/tools/psb_s_remote_mat.F90 @@ -98,20 +98,19 @@ Subroutine psb_ls_remote_mat(a,desc_a,b,info) integer(psb_ipk_) :: nnp, nrcvs, nsnds integer(psb_mpk_) :: icomm, minfo integer(psb_mpk_), allocatable :: brvindx(:), & - & rvsz(:), bsdindx(:),sdsz(:) + & rvsz(:), bsdindx(:),sdsz(:), sdsi(:), rvsi(:) 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(:) + integer(psb_ipk_), allocatable :: 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' + name='psb_s_remote_mat' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then info = psb_err_internal_error_ ; goto 9999 @@ -129,9 +128,9 @@ Subroutine psb_ls_remote_mat(a,desc_a,b,info) call b%free() - - Allocate(brvindx(np+1),& - & rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info) + + Allocate(rvsz(np),sdsz(np),sdsi(np),rvsi(np),brvindx(np+1),& + & bsdindx(np+1), acoo,stat=info) if (info /= psb_success_) then info=psb_err_alloc_dealloc_ @@ -140,71 +139,94 @@ Subroutine psb_ls_remote_mat(a,desc_a,b,info) 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) + 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' - 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 - + 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) - - - lnnz = max(iszr,iszs,lone) - lnc = a%get_ncols() - call acoo%allocate(lnr,lnc,lnnz) + icomm = desc_a%get_mpic() + sdsz(:)=0 + rvsz(:)=0 + sdsi(:)=0 + rvsi(:)=0 + ipx = 1 + brvindx(:) = 0 + bsdindx(:) = 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 + Do proc=0,np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + iszs = sum(sdsz) + iszr = sum(rvsz) + call acoo%allocate(desc_a%get_global_rows(),desc_a%get_global_cols(),iszr) + if (psb_errstatus_fatal()) then + write(0,*) 'Error from acoo%allocate ' + info = 4010 + goto 9999 + end if if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& & ' Send:',sdsz(:),' Receive:',rvsz(:) - + !write(debug_unit,*) me,' ',trim(name),': ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info) + !write(debug_unit,*) me,' ',trim(name),' iasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + !write(debug_unit,*) me,' ',trim(name),' jasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + !write(debug_unit,*) me,' ',trim(name),' 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 + do k=1, nz + proc = iprc(k) + sdsi(proc+1) = sdsi(proc+1) + 1 + !rvsi(proc) = rvsi(proc) + 1 + iasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ia(k) + jasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ja(k) + valsnd(bsdindx(proc+1)+sdsi(proc+1)) = a%val(k) + end do + do proc=0,np-1 + if (sdsi(proc+1) /= sdsz(proc+1)) & + & write(0,*) me,name,'Send mismacth ',sdsi(proc+1),sdsz(proc+1) + end do select case(psb_get_sp_a2av_alg()) case(psb_sp_a2av_smpl_triad_) @@ -218,7 +240,7 @@ Subroutine psb_ls_remote_mat(a,desc_a,b,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) & @@ -239,10 +261,9 @@ Subroutine psb_ls_remote_mat(a,desc_a,b,info) call psb_errpush(info,name,a_err='alltoallv') goto 9999 end if - + call acoo%set_nzeros(iszr) call acoo%mv_to_coo(b,info) - - + Deallocate(brvindx,bsdindx,rvsz,sdsz,& & iasnd,jasnd,valsnd,stat=info) if (debug_level >= psb_debug_outer_)& diff --git a/base/tools/psb_sspalloc.f90 b/base/tools/psb_sspalloc.f90 index 2d79b0ba..14e784da 100644 --- a/base/tools/psb_sspalloc.f90 +++ b/base/tools/psb_sspalloc.f90 @@ -114,16 +114,16 @@ subroutine psb_sspalloc(a, desc_a, info, nnz, bldmode) goto 9999 end if - write(0,*) name,'Setting a%remote_build ',& - & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ +!!$ write(0,*) name,'Setting a%remote_build ',& +!!$ & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ a%remote_build = bldmode_ select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed - write(0,*) name,' matbld_noremote_ nothing needed' + !write(0,*) name,' matbld_noremote_ nothing needed' case (psb_matbld_remote_) - write(0,*) name,' matbld_remote_ start ' + !write(0,*) name,' matbld_remote_ start ' allocate(a%rmta) nnzrmt_ = max(100,(nnz_/100)) call a%rmta%allocate(m,n,nnzrmt_) diff --git a/base/tools/psb_sspasb.f90 b/base/tools/psb_sspasb.f90 index fa9a2a49..9249c59c 100644 --- a/base/tools/psb_sspasb.f90 +++ b/base/tools/psb_sspasb.f90 @@ -106,26 +106,38 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) ! First case: we come from a fresh build. ! - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - call a%set_nrows(n_row) - call a%set_ncols(n_col) - end if - - if (a%is_bld()) then - select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) me,name,' 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 - + integer(psb_ipk_), allocatable :: ila(:), jla(:) + integer(psb_ipk_) :: nz, nzt,k call psb_remote_mat(a%rmta,desc_a,a_add,info) - - end block + nz = a_add%get_nzeros() +!!$ write(0,*) me,name,' Nz to be added',nz + nzt = nz + call psb_sum(ctxt,nzt) + if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nz > 0) then + ! + ! Should we check for new indices here? + ! + call psb_realloc(nz,ila,info) + call psb_realloc(nz,jla,info) + call desc_a%indxmap%g2l(a_add%ia(1:nz),ila(1:nz),info,owned=.true.) + if (info == 0) call desc_a%indxmap%g2l_ins(a_add%ja(1:nz),jla(1:nz),info) + !write(0,*) me,name,' Check before insert',a%get_nzeros() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) + !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz + end if + if (nzt > 0) call psb_cdasb(desc_a,info) + end block end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_sspins.F90 b/base/tools/psb_sspins.F90 index c4abb9f6..e49e7423 100644 --- a/base/tools/psb_sspins.F90 +++ b/base/tools/psb_sspins.F90 @@ -164,9 +164,9 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) do i=1,nz if (ila(i)<0) then k=k+1 - lila(k) = ia(k) - ljla(k) = ja(k) - lval(k) = val(k) + lila(k) = ia(i) + ljla(k) = ja(i) + lval(k) = val(i) end if end do if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl diff --git a/base/tools/psb_z_remote_mat.F90 b/base/tools/psb_z_remote_mat.F90 index 1acc57b2..7b2ade7c 100644 --- a/base/tools/psb_z_remote_mat.F90 +++ b/base/tools/psb_z_remote_mat.F90 @@ -98,20 +98,19 @@ Subroutine psb_lz_remote_mat(a,desc_a,b,info) integer(psb_ipk_) :: nnp, nrcvs, nsnds integer(psb_mpk_) :: icomm, minfo integer(psb_mpk_), allocatable :: brvindx(:), & - & rvsz(:), bsdindx(:),sdsz(:) + & rvsz(:), bsdindx(:),sdsz(:), sdsi(:), rvsi(:) 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(:) + integer(psb_ipk_), allocatable :: 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' + name='psb_z_remote_mat' call psb_erractionsave(err_act) if (psb_errstatus_fatal()) then info = psb_err_internal_error_ ; goto 9999 @@ -129,9 +128,9 @@ Subroutine psb_lz_remote_mat(a,desc_a,b,info) call b%free() - - Allocate(brvindx(np+1),& - & rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info) + + Allocate(rvsz(np),sdsz(np),sdsi(np),rvsi(np),brvindx(np+1),& + & bsdindx(np+1), acoo,stat=info) if (info /= psb_success_) then info=psb_err_alloc_dealloc_ @@ -140,71 +139,94 @@ Subroutine psb_lz_remote_mat(a,desc_a,b,info) 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) + 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' - 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 - + 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) - - - lnnz = max(iszr,iszs,lone) - lnc = a%get_ncols() - call acoo%allocate(lnr,lnc,lnnz) + icomm = desc_a%get_mpic() + sdsz(:)=0 + rvsz(:)=0 + sdsi(:)=0 + rvsi(:)=0 + ipx = 1 + brvindx(:) = 0 + bsdindx(:) = 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 + Do proc=0,np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + iszs = sum(sdsz) + iszr = sum(rvsz) + call acoo%allocate(desc_a%get_global_rows(),desc_a%get_global_cols(),iszr) + if (psb_errstatus_fatal()) then + write(0,*) 'Error from acoo%allocate ' + info = 4010 + goto 9999 + end if if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& & ' Send:',sdsz(:),' Receive:',rvsz(:) - + !write(debug_unit,*) me,' ',trim(name),': ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info) + !write(debug_unit,*) me,' ',trim(name),' iasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + !write(debug_unit,*) me,' ',trim(name),' jasnd: ',info if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + !write(debug_unit,*) me,' ',trim(name),' 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 + do k=1, nz + proc = iprc(k) + sdsi(proc+1) = sdsi(proc+1) + 1 + !rvsi(proc) = rvsi(proc) + 1 + iasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ia(k) + jasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ja(k) + valsnd(bsdindx(proc+1)+sdsi(proc+1)) = a%val(k) + end do + do proc=0,np-1 + if (sdsi(proc+1) /= sdsz(proc+1)) & + & write(0,*) me,name,'Send mismacth ',sdsi(proc+1),sdsz(proc+1) + end do select case(psb_get_sp_a2av_alg()) case(psb_sp_a2av_smpl_triad_) @@ -218,7 +240,7 @@ Subroutine psb_lz_remote_mat(a,desc_a,b,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) & @@ -239,10 +261,9 @@ Subroutine psb_lz_remote_mat(a,desc_a,b,info) call psb_errpush(info,name,a_err='alltoallv') goto 9999 end if - + call acoo%set_nzeros(iszr) call acoo%mv_to_coo(b,info) - - + Deallocate(brvindx,bsdindx,rvsz,sdsz,& & iasnd,jasnd,valsnd,stat=info) if (debug_level >= psb_debug_outer_)& diff --git a/base/tools/psb_zspalloc.f90 b/base/tools/psb_zspalloc.f90 index 1d76d66c..3b0fc7d2 100644 --- a/base/tools/psb_zspalloc.f90 +++ b/base/tools/psb_zspalloc.f90 @@ -114,16 +114,16 @@ subroutine psb_zspalloc(a, desc_a, info, nnz, bldmode) goto 9999 end if - write(0,*) name,'Setting a%remote_build ',& - & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ +!!$ write(0,*) name,'Setting a%remote_build ',& +!!$ & bldmode_,psb_matbld_noremote_,psb_matbld_remote_ a%remote_build = bldmode_ select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed - write(0,*) name,' matbld_noremote_ nothing needed' + !write(0,*) name,' matbld_noremote_ nothing needed' case (psb_matbld_remote_) - write(0,*) name,' matbld_remote_ start ' + !write(0,*) name,' matbld_remote_ start ' allocate(a%rmta) nnzrmt_ = max(100,(nnz_/100)) call a%rmta%allocate(m,n,nnzrmt_) diff --git a/base/tools/psb_zspasb.f90 b/base/tools/psb_zspasb.f90 index 67353739..d92daa2c 100644 --- a/base/tools/psb_zspasb.f90 +++ b/base/tools/psb_zspasb.f90 @@ -106,26 +106,38 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) ! First case: we come from a fresh build. ! - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() - call a%set_nrows(n_row) - call a%set_ncols(n_col) - end if - - if (a%is_bld()) then - select case(a%remote_build) case (psb_matbld_noremote_) ! nothing needed case (psb_matbld_remote_) - write(0,*) me,name,' 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 - + integer(psb_ipk_), allocatable :: ila(:), jla(:) + integer(psb_ipk_) :: nz, nzt,k call psb_remote_mat(a%rmta,desc_a,a_add,info) - - end block + nz = a_add%get_nzeros() +!!$ write(0,*) me,name,' Nz to be added',nz + nzt = nz + call psb_sum(ctxt,nzt) + if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nz > 0) then + ! + ! Should we check for new indices here? + ! + call psb_realloc(nz,ila,info) + call psb_realloc(nz,jla,info) + call desc_a%indxmap%g2l(a_add%ia(1:nz),ila(1:nz),info,owned=.true.) + if (info == 0) call desc_a%indxmap%g2l_ins(a_add%ja(1:nz),jla(1:nz),info) + !write(0,*) me,name,' Check before insert',a%get_nzeros() + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) + !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz + end if + if (nzt > 0) call psb_cdasb(desc_a,info) + end block end select call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) diff --git a/base/tools/psb_zspins.F90 b/base/tools/psb_zspins.F90 index 2a654e14..835a4d04 100644 --- a/base/tools/psb_zspins.F90 +++ b/base/tools/psb_zspins.F90 @@ -164,9 +164,9 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local) do i=1,nz if (ila(i)<0) then k=k+1 - lila(k) = ia(k) - ljla(k) = ja(k) - lval(k) = val(k) + lila(k) = ia(i) + ljla(k) = ja(i) + lval(k) = val(i) end if end do if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl diff --git a/test/pargen/psb_d_pde3d.F90 b/test/pargen/psb_d_pde3d.F90 index 483e6d9c..2384d410 100644 --- a/test/pargen/psb_d_pde3d.F90 +++ b/test/pargen/psb_d_pde3d.F90 @@ -679,9 +679,9 @@ contains t1 = psb_wtime() if (info == psb_success_) then if (present(amold)) then - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + call psb_spasb(a,desc_a,info,dupl=psb_dupl_add_,mold=amold) else - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + call psb_spasb(a,desc_a,info,dupl=psb_dupl_add_,afmt=afmt) end if end if call psb_barrier(ctxt) diff --git a/test/pargen/psb_s_pde3d.F90 b/test/pargen/psb_s_pde3d.F90 index b27f8818..96dbd5db 100644 --- a/test/pargen/psb_s_pde3d.F90 +++ b/test/pargen/psb_s_pde3d.F90 @@ -679,9 +679,9 @@ contains t1 = psb_wtime() if (info == psb_success_) then if (present(amold)) then - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) + call psb_spasb(a,desc_a,info,dupl=psb_dupl_add_,mold=amold) else - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) + call psb_spasb(a,desc_a,info,dupl=psb_dupl_add_,afmt=afmt) end if end if call psb_barrier(ctxt) From 7d150e2eca2eea0cfabea3757eb3cfa4f7b453d8 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sun, 6 Feb 2022 18:34:43 +0100 Subject: [PATCH 05/28] Fix assembly procedure for remote build --- base/tools/psb_cspasb.f90 | 1 + base/tools/psb_dspasb.f90 | 1 + base/tools/psb_sspasb.f90 | 1 + base/tools/psb_zspasb.f90 | 1 + 4 files changed, 4 insertions(+) diff --git a/base/tools/psb_cspasb.f90 b/base/tools/psb_cspasb.f90 index 99c8990e..6cec5fa7 100644 --- a/base/tools/psb_cspasb.f90 +++ b/base/tools/psb_cspasb.f90 @@ -139,6 +139,7 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) end block end select + call a%set_ncols(desc_a%get_local_cols()) call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index c9d3cd08..422a8661 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -139,6 +139,7 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) end block end select + call a%set_ncols(desc_a%get_local_cols()) call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then diff --git a/base/tools/psb_sspasb.f90 b/base/tools/psb_sspasb.f90 index 9249c59c..d6de2136 100644 --- a/base/tools/psb_sspasb.f90 +++ b/base/tools/psb_sspasb.f90 @@ -139,6 +139,7 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) end block end select + call a%set_ncols(desc_a%get_local_cols()) call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then diff --git a/base/tools/psb_zspasb.f90 b/base/tools/psb_zspasb.f90 index d92daa2c..b343426d 100644 --- a/base/tools/psb_zspasb.f90 +++ b/base/tools/psb_zspasb.f90 @@ -139,6 +139,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) end block end select + call a%set_ncols(desc_a%get_local_cols()) call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) else if (a%is_upd()) then From 4a7f9d786d60746547e2249bd1f8409b6b818ac8 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 7 Feb 2022 09:10:30 +0100 Subject: [PATCH 06/28] Keep track of inner integer class for descriptor reassembly --- base/tools/psb_cspasb.f90 | 9 +++++++-- base/tools/psb_dspasb.f90 | 9 +++++++-- base/tools/psb_sspasb.f90 | 9 +++++++-- base/tools/psb_zspasb.f90 | 9 +++++++-- 4 files changed, 28 insertions(+), 8 deletions(-) diff --git a/base/tools/psb_cspasb.f90 b/base/tools/psb_cspasb.f90 index 6cec5fa7..390bdcaf 100644 --- a/base/tools/psb_cspasb.f90 +++ b/base/tools/psb_cspasb.f90 @@ -68,6 +68,7 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err + class(psb_i_base_vect_type), allocatable :: ivm info = psb_success_ name = 'psb_spasb' @@ -120,7 +121,10 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) !!$ write(0,*) me,name,' Nz to be added',nz nzt = nz call psb_sum(ctxt,nzt) - if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nzt>0) then + allocate(ivm, mold=desc_a%v_halo_index%v) + call psb_cd_reinit(desc_a, info) + end if if (nz > 0) then ! ! Should we check for new indices here? @@ -132,10 +136,11 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) !write(0,*) me,name,' Check before insert',a%get_nzeros() n_row = desc_a%get_local_rows() n_col = desc_a%get_local_cols() + call a%set_ncols(desc_a%get_local_cols()) call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz end if - if (nzt > 0) call psb_cdasb(desc_a,info) + if (nzt > 0) call psb_cdasb(desc_a,info,mold=ivm) end block end select diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index 422a8661..d82f5efe 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -68,6 +68,7 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err + class(psb_i_base_vect_type), allocatable :: ivm info = psb_success_ name = 'psb_spasb' @@ -120,7 +121,10 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) !!$ write(0,*) me,name,' Nz to be added',nz nzt = nz call psb_sum(ctxt,nzt) - if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nzt>0) then + allocate(ivm, mold=desc_a%v_halo_index%v) + call psb_cd_reinit(desc_a, info) + end if if (nz > 0) then ! ! Should we check for new indices here? @@ -132,10 +136,11 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) !write(0,*) me,name,' Check before insert',a%get_nzeros() n_row = desc_a%get_local_rows() n_col = desc_a%get_local_cols() + call a%set_ncols(desc_a%get_local_cols()) call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz end if - if (nzt > 0) call psb_cdasb(desc_a,info) + if (nzt > 0) call psb_cdasb(desc_a,info,mold=ivm) end block end select diff --git a/base/tools/psb_sspasb.f90 b/base/tools/psb_sspasb.f90 index d6de2136..87372d01 100644 --- a/base/tools/psb_sspasb.f90 +++ b/base/tools/psb_sspasb.f90 @@ -68,6 +68,7 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err + class(psb_i_base_vect_type), allocatable :: ivm info = psb_success_ name = 'psb_spasb' @@ -120,7 +121,10 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) !!$ write(0,*) me,name,' Nz to be added',nz nzt = nz call psb_sum(ctxt,nzt) - if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nzt>0) then + allocate(ivm, mold=desc_a%v_halo_index%v) + call psb_cd_reinit(desc_a, info) + end if if (nz > 0) then ! ! Should we check for new indices here? @@ -132,10 +136,11 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) !write(0,*) me,name,' Check before insert',a%get_nzeros() n_row = desc_a%get_local_rows() n_col = desc_a%get_local_cols() + call a%set_ncols(desc_a%get_local_cols()) call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz end if - if (nzt > 0) call psb_cdasb(desc_a,info) + if (nzt > 0) call psb_cdasb(desc_a,info,mold=ivm) end block end select diff --git a/base/tools/psb_zspasb.f90 b/base/tools/psb_zspasb.f90 index b343426d..9c42ae85 100644 --- a/base/tools/psb_zspasb.f90 +++ b/base/tools/psb_zspasb.f90 @@ -68,6 +68,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err + class(psb_i_base_vect_type), allocatable :: ivm info = psb_success_ name = 'psb_spasb' @@ -120,7 +121,10 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) !!$ write(0,*) me,name,' Nz to be added',nz nzt = nz call psb_sum(ctxt,nzt) - if (nzt>0) call psb_cd_reinit(desc_a, info) + if (nzt>0) then + allocate(ivm, mold=desc_a%v_halo_index%v) + call psb_cd_reinit(desc_a, info) + end if if (nz > 0) then ! ! Should we check for new indices here? @@ -132,10 +136,11 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) !write(0,*) me,name,' Check before insert',a%get_nzeros() n_row = desc_a%get_local_rows() n_col = desc_a%get_local_cols() + call a%set_ncols(desc_a%get_local_cols()) call a%csput(nz,ila,jla,a_add%val,ione,n_row,ione,n_col,info) !write(0,*) me,name,' Check after insert',a%get_nzeros(),nz end if - if (nzt > 0) call psb_cdasb(desc_a,info) + if (nzt > 0) call psb_cdasb(desc_a,info,mold=ivm) end block end select From fef69bf799088d1ad42e1eac44b845ebc3f735ef Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 7 Feb 2022 16:54:37 +0100 Subject: [PATCH 07/28] Update docs for remote build. --- docs/html/index.html | 2 +- docs/html/userhtml.html | 2 +- docs/html/userhtmlsu36.html | 35 +- docs/html/userhtmlsu37.html | 14 +- docs/html/userhtmlsu38.html | 27 +- docs/html/userhtmlsu39.html | 12 +- docs/html/userhtmlsu40.html | 14 +- docs/html/userhtmlsu41.html | 12 +- docs/html/userhtmlsu42.html | 14 +- docs/html/userhtmlsu43.html | 12 +- docs/html/userhtmlsu44.html | 12 +- docs/html/userhtmlsu45.html | 12 +- docs/html/userhtmlsu46.html | 14 +- docs/html/userhtmlsu47.html | 12 +- docs/html/userhtmlsu48.html | 14 +- docs/html/userhtmlsu49.html | 14 +- docs/html/userhtmlsu50.html | 14 +- docs/html/userhtmlsu51.html | 14 +- docs/html/userhtmlsu52.html | 14 +- docs/html/userhtmlsu53.html | 14 +- docs/html/userhtmlsu54.html | 14 +- docs/html/userhtmlsu55.html | 12 +- docs/html/userhtmlsu56.html | 16 +- docs/psblas-3.7.pdf | 3101 ++++++++++++++++++----------------- docs/src/toolsrout.tex | 19 +- docs/src/userguide.tex | 4 +- docs/src/userhtml.tex | 2 +- 27 files changed, 1753 insertions(+), 1693 deletions(-) diff --git a/docs/html/index.html b/docs/html/index.html index 2906d9ce..ece52ac5 100644 --- a/docs/html/index.html +++ b/docs/html/index.html @@ -21,7 +21,7 @@ class="cmbx-10">Salvatore Filippone
Alfredo Buttari
Software version: 3.7.0.1
May 11th, 2021 +class="newline" />Feb 14th, 2022 diff --git a/docs/html/userhtml.html b/docs/html/userhtml.html index 2906d9ce..ece52ac5 100644 --- a/docs/html/userhtml.html +++ b/docs/html/userhtml.html @@ -21,7 +21,7 @@ class="cmbx-10">Salvatore Filippone
Alfredo Buttari
Software version: 3.7.0.1
May 11th, 2021 +class="newline" />Feb 14th, 2022 diff --git a/docs/html/userhtmlsu36.html b/docs/html/userhtmlsu36.html index 9616adcb..dfb56550 100644 --- a/docs/html/userhtmlsu36.html +++ b/docs/html/userhtmlsu36.html @@ -22,7 +22,7 @@ href="userhtmlsu32.html#userhtmlsu39.html" >up]

-call psb_spall(a, desc_a, info, nnz)
+call psb_spall(a, desc_a, info, nnz, buildmode)
 

@@ -61,8 +61,25 @@ class="newline" />Type: optional.
Intent: in.
Specified as: an integer value. -

+class="newline" />Specified as: an integer value. +

+buildmode
Whether to keep track of matrix entries that do not belong to the + current process.
Scope: global.
Type: optional.
Intent: in.
Specified as: + an integer value psb_matbld_noremote_, psb_matbld_remote_. Default: + psb_matbld_noremote_.
+

On Return
psb_Tspmat_type. + + +
info
required
Intent: out.
An integer value; 0 means no error has been detected.
-

Notes

  1. On exit from this routine the sparse matrix is in the build state. - - -
  2. The descriptor may be in either the build or assembled state. @@ -117,12 +134,12 @@ class="cmmi-10">nnz in the -