From b534a24b7a40232216473b2b8b261aeecb782a60 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 14 Oct 2020 17:32:06 +0200 Subject: [PATCH] New ptap_bld function. Modify ptap not to rebuild Restr. --- mlprec/amg_c_base_aggregator_mod.f90 | 21 +- mlprec/amg_d_base_aggregator_mod.f90 | 21 +- mlprec/amg_s_base_aggregator_mod.f90 | 21 +- mlprec/amg_z_base_aggregator_mod.f90 | 21 +- mlprec/impl/aggregator/Makefile | 8 +- mlprec/impl/aggregator/amg_c_ptap.f90 | 571 +------------- mlprec/impl/aggregator/amg_c_ptap_bld.f90 | 697 ++++++++++++++++++ .../aggregator/amg_caggrmat_nosmth_bld.f90 | 2 +- .../impl/aggregator/amg_caggrmat_smth_bld.f90 | 2 +- mlprec/impl/aggregator/amg_d_ptap.f90 | 571 +------------- mlprec/impl/aggregator/amg_d_ptap_bld.f90 | 697 ++++++++++++++++++ .../aggregator/amg_daggrmat_nosmth_bld.f90 | 2 +- .../impl/aggregator/amg_daggrmat_smth_bld.f90 | 2 +- mlprec/impl/aggregator/amg_s_ptap.f90 | 571 +------------- mlprec/impl/aggregator/amg_s_ptap_bld.f90 | 697 ++++++++++++++++++ .../aggregator/amg_saggrmat_nosmth_bld.f90 | 2 +- .../impl/aggregator/amg_saggrmat_smth_bld.f90 | 2 +- mlprec/impl/aggregator/amg_z_ptap.f90 | 571 +------------- mlprec/impl/aggregator/amg_z_ptap_bld.f90 | 697 ++++++++++++++++++ .../aggregator/amg_zaggrmat_nosmth_bld.f90 | 2 +- .../impl/aggregator/amg_zaggrmat_smth_bld.f90 | 2 +- 21 files changed, 2964 insertions(+), 2216 deletions(-) create mode 100644 mlprec/impl/aggregator/amg_c_ptap_bld.f90 create mode 100644 mlprec/impl/aggregator/amg_d_ptap_bld.f90 create mode 100644 mlprec/impl/aggregator/amg_s_ptap_bld.f90 create mode 100644 mlprec/impl/aggregator/amg_z_ptap_bld.f90 diff --git a/mlprec/amg_c_base_aggregator_mod.f90 b/mlprec/amg_c_base_aggregator_mod.f90 index ec0b2915..3d6485f4 100644 --- a/mlprec/amg_c_base_aggregator_mod.f90 +++ b/mlprec/amg_c_base_aggregator_mod.f90 @@ -119,8 +119,8 @@ module amg_c_base_aggregator_mod end subroutine amg_c_soc_map_bld end interface - interface amg_ptap - subroutine amg_c_ptap(a_csr,desc_a,nlaggr,parms,ac,& + interface amg_ptap_bld + subroutine amg_c_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_cprol,coo_restr,info,desc_ax) import :: psb_c_csr_sparse_mat, psb_cspmat_type, psb_desc_type, & & psb_c_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ @@ -134,6 +134,23 @@ module amg_c_base_aggregator_mod type(psb_cspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info type(psb_desc_type), intent(inout), optional :: desc_ax + end subroutine amg_c_ptap_bld + end interface amg_ptap_bld + + interface amg_ptap + subroutine amg_c_ptap(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_cprol,coo_restr,info) + import :: psb_c_csr_sparse_mat, psb_cspmat_type, psb_desc_type, & + & psb_c_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ + implicit none + type(psb_c_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(in) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_c_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_cprol + type(psb_cspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info end subroutine amg_c_ptap end interface amg_ptap diff --git a/mlprec/amg_d_base_aggregator_mod.f90 b/mlprec/amg_d_base_aggregator_mod.f90 index 8dcf8b65..5ec0205c 100644 --- a/mlprec/amg_d_base_aggregator_mod.f90 +++ b/mlprec/amg_d_base_aggregator_mod.f90 @@ -119,8 +119,8 @@ module amg_d_base_aggregator_mod end subroutine amg_d_soc_map_bld end interface - interface amg_ptap - subroutine amg_d_ptap(a_csr,desc_a,nlaggr,parms,ac,& + interface amg_ptap_bld + subroutine amg_d_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_cprol,coo_restr,info,desc_ax) import :: psb_d_csr_sparse_mat, psb_dspmat_type, psb_desc_type, & & psb_d_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ @@ -134,6 +134,23 @@ module amg_d_base_aggregator_mod type(psb_dspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info type(psb_desc_type), intent(inout), optional :: desc_ax + end subroutine amg_d_ptap_bld + end interface amg_ptap_bld + + interface amg_ptap + subroutine amg_d_ptap(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_cprol,coo_restr,info) + import :: psb_d_csr_sparse_mat, psb_dspmat_type, psb_desc_type, & + & psb_d_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ + implicit none + type(psb_d_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(in) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_d_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_cprol + type(psb_dspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info end subroutine amg_d_ptap end interface amg_ptap diff --git a/mlprec/amg_s_base_aggregator_mod.f90 b/mlprec/amg_s_base_aggregator_mod.f90 index 97bf7dd2..a3545bdf 100644 --- a/mlprec/amg_s_base_aggregator_mod.f90 +++ b/mlprec/amg_s_base_aggregator_mod.f90 @@ -119,8 +119,8 @@ module amg_s_base_aggregator_mod end subroutine amg_s_soc_map_bld end interface - interface amg_ptap - subroutine amg_s_ptap(a_csr,desc_a,nlaggr,parms,ac,& + interface amg_ptap_bld + subroutine amg_s_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_cprol,coo_restr,info,desc_ax) import :: psb_s_csr_sparse_mat, psb_sspmat_type, psb_desc_type, & & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ @@ -134,6 +134,23 @@ module amg_s_base_aggregator_mod type(psb_sspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info type(psb_desc_type), intent(inout), optional :: desc_ax + end subroutine amg_s_ptap_bld + end interface amg_ptap_bld + + interface amg_ptap + subroutine amg_s_ptap(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_cprol,coo_restr,info) + import :: psb_s_csr_sparse_mat, psb_sspmat_type, psb_desc_type, & + & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ + implicit none + type(psb_s_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(in) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_cprol + type(psb_sspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info end subroutine amg_s_ptap end interface amg_ptap diff --git a/mlprec/amg_z_base_aggregator_mod.f90 b/mlprec/amg_z_base_aggregator_mod.f90 index 9bc7068c..5d039ef6 100644 --- a/mlprec/amg_z_base_aggregator_mod.f90 +++ b/mlprec/amg_z_base_aggregator_mod.f90 @@ -119,8 +119,8 @@ module amg_z_base_aggregator_mod end subroutine amg_z_soc_map_bld end interface - interface amg_ptap - subroutine amg_z_ptap(a_csr,desc_a,nlaggr,parms,ac,& + interface amg_ptap_bld + subroutine amg_z_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_cprol,coo_restr,info,desc_ax) import :: psb_z_csr_sparse_mat, psb_zspmat_type, psb_desc_type, & & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ @@ -134,6 +134,23 @@ module amg_z_base_aggregator_mod type(psb_zspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info type(psb_desc_type), intent(inout), optional :: desc_ax + end subroutine amg_z_ptap_bld + end interface amg_ptap_bld + + interface amg_ptap + subroutine amg_z_ptap(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_cprol,coo_restr,info) + import :: psb_z_csr_sparse_mat, psb_zspmat_type, psb_desc_type, & + & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ + implicit none + type(psb_z_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(in) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_cprol + type(psb_zspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info end subroutine amg_z_ptap end interface amg_ptap diff --git a/mlprec/impl/aggregator/Makefile b/mlprec/impl/aggregator/Makefile index 71ebe6f3..a7dce5be 100644 --- a/mlprec/impl/aggregator/Makefile +++ b/mlprec/impl/aggregator/Makefile @@ -14,7 +14,7 @@ amg_s_dec_aggregator_mat_bld.o \ amg_s_dec_aggregator_tprol.o \ amg_s_symdec_aggregator_tprol.o \ amg_s_map_to_tprol.o amg_s_soc1_map_bld.o amg_s_soc2_map_bld.o\ -amg_s_ptap.o \ +amg_s_ptap.o amg_s_ptap_bld.o \ amg_saggrmat_minnrg_bld.o\ amg_saggrmat_nosmth_bld.o amg_saggrmat_smth_bld.o \ amg_d_dec_aggregator_mat_asb.o \ @@ -22,7 +22,7 @@ amg_d_dec_aggregator_mat_bld.o \ amg_d_dec_aggregator_tprol.o \ amg_d_symdec_aggregator_tprol.o \ amg_d_map_to_tprol.o amg_d_soc1_map_bld.o amg_d_soc2_map_bld.o \ -amg_d_ptap.o \ +amg_d_ptap.o amg_d_ptap_bld.o \ amg_daggrmat_minnrg_bld.o \ amg_daggrmat_nosmth_bld.o amg_daggrmat_smth_bld.o \ amg_c_dec_aggregator_mat_asb.o \ @@ -30,7 +30,7 @@ amg_c_dec_aggregator_mat_bld.o \ amg_c_dec_aggregator_tprol.o \ amg_c_symdec_aggregator_tprol.o \ amg_c_map_to_tprol.o amg_c_soc1_map_bld.o amg_c_soc2_map_bld.o\ -amg_c_ptap.o \ +amg_c_ptap.o amg_c_ptap_bld.o \ amg_caggrmat_minnrg_bld.o\ amg_caggrmat_nosmth_bld.o amg_caggrmat_smth_bld.o \ amg_z_dec_aggregator_mat_asb.o \ @@ -38,7 +38,7 @@ amg_z_dec_aggregator_mat_bld.o \ amg_z_dec_aggregator_tprol.o \ amg_z_symdec_aggregator_tprol.o \ amg_z_map_to_tprol.o amg_z_soc1_map_bld.o amg_z_soc2_map_bld.o\ -amg_z_ptap.o \ +amg_z_ptap.o amg_z_ptap_bld.o \ amg_zaggrmat_minnrg_bld.o\ amg_zaggrmat_nosmth_bld.o amg_zaggrmat_smth_bld.o diff --git a/mlprec/impl/aggregator/amg_c_ptap.f90 b/mlprec/impl/aggregator/amg_c_ptap.f90 index 035a6b89..42ed4202 100644 --- a/mlprec/impl/aggregator/amg_c_ptap.f90 +++ b/mlprec/impl/aggregator/amg_c_ptap.f90 @@ -32,11 +32,14 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -! File: amg_daggrmat_nosmth_bld.F90 +! File: amg_ptap.F90 ! +! This routine computes the triple product : +! AC = P^T A P +! P^T is stored in Restr ! subroutine amg_c_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) + & coo_prol,desc_ac,coo_restr,info) use psb_base_mod use amg_c_inner_mod use amg_c_base_aggregator_mod, amg_protect_name => amg_c_ptap @@ -51,7 +54,6 @@ subroutine amg_c_ptap(a_csr,desc_a,nlaggr,parms,ac,& type(psb_desc_type), intent(inout) :: desc_ac type(psb_cspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax ! Local variables integer(psb_ipk_) :: err_act @@ -64,8 +66,7 @@ subroutine amg_c_ptap(a_csr,desc_a,nlaggr,parms,ac,& integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 + logical, parameter :: debug=.false. name='amg_ptap' if(psb_get_errstatus().ne.0) return @@ -82,9 +83,6 @@ subroutine amg_c_ptap(a_csr,desc_a,nlaggr,parms,ac,& nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) @@ -106,542 +104,38 @@ subroutine amg_c_ptap(a_csr,desc_a,nlaggr,parms,ac,& & desc_ac%get_local_rows(),desc_a%get_local_cols() if (debug) flush(0) - if (do_timings) call psb_tic(idx_spspmm) call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) if (debug) write(0,*) me,trim(name),' Done AxPROL ',& & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& & desc_ac%get_local_rows(),desc_ac%get_local_cols() - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - call coo_prol%cp_to_coo(coo_restr,info) - call coo_restr%set_ncols(desc_ac%get_local_cols()) - call coo_restr%set_nrows(desc_a%get_local_rows()) - call psb_c_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_c_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_c_ptap - -subroutine amg_c_lc_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_c_inner_mod - use amg_c_base_aggregator_mod !, amg_protect_name => amg_c_lc_ptap - implicit none - - ! Arguments - type(psb_c_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_sml_parms), intent(inout) :: parms - type(psb_lc_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_lcspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_c_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_ifmt(csr_prol,info) + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - type(psb_c_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_c_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_lcoo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_lcoo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_lc_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_c_lc_ptap - -subroutine amg_lc_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_c_inner_mod - use amg_c_base_aggregator_mod!, amg_protect_name => amg_lc_ptap - implicit none - - ! Arguments - type(psb_lc_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_sml_parms), intent(inout) :: parms - type(psb_lc_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_lcspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_lc_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr - - ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_fmt(csr_prol,info) - - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. + & 'starting sphalo/ rwxtd' - if (present(desc_ax)) then - block - type(psb_c_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_c_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) - write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& - & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() - if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) call ac_csr%set_nrows(desc_ac%get_local_rows()) call ac_csr%set_ncols(desc_ac%get_local_cols()) @@ -653,12 +147,8 @@ subroutine amg_lc_ptap(a_csr,desc_a,nlaggr,parms,ac,& ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() call coo_prol%set_ncols(desc_ac%get_local_cols()) - !call coo_restr%mv_from_ifmt(csr_restr,info) -!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) -!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done ptap ' @@ -669,21 +159,4 @@ subroutine amg_lc_ptap(a_csr,desc_a,nlaggr,parms,ac,& 9999 call psb_error_handler(err_act) return - -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_lc_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_lc_ptap +end subroutine amg_c_ptap diff --git a/mlprec/impl/aggregator/amg_c_ptap_bld.f90 b/mlprec/impl/aggregator/amg_c_ptap_bld.f90 new file mode 100644 index 00000000..fa1fa19a --- /dev/null +++ b/mlprec/impl/aggregator/amg_c_ptap_bld.f90 @@ -0,0 +1,697 @@ +! +! +! MLD2P4 Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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: amg_ptap_bld.F90 +! +! This routine does two things: +! 1. Computes the transpose of the prolongator Prol and stores +! it explicitly into Restr +! 2. Computes AC = P^T A P +! +! Yes, doing two things at once is against received wisdom, but we are +! trying to save time in the build procedure; maybe we will revisit this +! decision. +! +subroutine amg_c_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_c_inner_mod + use amg_c_base_aggregator_mod, amg_protect_name => amg_c_ptap_bld + implicit none + + ! Arguments + type(psb_c_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_c_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_cspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_c_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 + integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + call coo_prol%cp_to_coo(coo_restr,info) + call coo_restr%set_ncols(desc_ac%get_local_cols()) + call coo_restr%set_nrows(desc_a%get_local_rows()) + call psb_c_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_c_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_c_ptap_bld + +subroutine amg_c_lc_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_c_inner_mod + use amg_c_base_aggregator_mod !, amg_protect_name => amg_c_lc_ptap_bld + implicit none + + ! Arguments + type(psb_c_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_lc_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_lcspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_c_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_ifmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_c_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_c_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_lcoo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_lcoo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_lc_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_c_lc_ptap_bld + +subroutine amg_lc_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_c_inner_mod + use amg_c_base_aggregator_mod!, amg_protect_name => amg_lc_ptap_bld + implicit none + + ! Arguments + type(psb_lc_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_lc_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_lcspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_lc_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_c_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_c_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) + write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& + & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() + if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + !call coo_restr%mv_from_ifmt(csr_restr,info) +!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) +!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_lc_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_lc_ptap_bld diff --git a/mlprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 b/mlprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 index def8bb08..af00bf1c 100644 --- a/mlprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 +++ b/mlprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 @@ -160,7 +160,7 @@ subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& call psb_cdasb(desc_ac,info) call psb_cd_reinit(desc_ac,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call coo_restr%set_nrows(desc_ac%get_local_rows()) diff --git a/mlprec/impl/aggregator/amg_caggrmat_smth_bld.f90 b/mlprec/impl/aggregator/amg_caggrmat_smth_bld.f90 index 53fee0f2..b1011ea1 100644 --- a/mlprec/impl/aggregator/amg_caggrmat_smth_bld.f90 +++ b/mlprec/impl/aggregator/amg_caggrmat_smth_bld.f90 @@ -286,7 +286,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& nzl = acsr1%get_nzeros() call acsr1%mv_to_coo(coo_prol,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call op_prol%mv_from(coo_prol) diff --git a/mlprec/impl/aggregator/amg_d_ptap.f90 b/mlprec/impl/aggregator/amg_d_ptap.f90 index 9f3c8c54..fbb7ac65 100644 --- a/mlprec/impl/aggregator/amg_d_ptap.f90 +++ b/mlprec/impl/aggregator/amg_d_ptap.f90 @@ -32,11 +32,14 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -! File: amg_daggrmat_nosmth_bld.F90 +! File: amg_ptap.F90 ! +! This routine computes the triple product : +! AC = P^T A P +! P^T is stored in Restr ! subroutine amg_d_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) + & coo_prol,desc_ac,coo_restr,info) use psb_base_mod use amg_d_inner_mod use amg_d_base_aggregator_mod, amg_protect_name => amg_d_ptap @@ -51,7 +54,6 @@ subroutine amg_d_ptap(a_csr,desc_a,nlaggr,parms,ac,& type(psb_desc_type), intent(inout) :: desc_ac type(psb_dspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax ! Local variables integer(psb_ipk_) :: err_act @@ -64,8 +66,7 @@ subroutine amg_d_ptap(a_csr,desc_a,nlaggr,parms,ac,& integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 + logical, parameter :: debug=.false. name='amg_ptap' if(psb_get_errstatus().ne.0) return @@ -82,9 +83,6 @@ subroutine amg_d_ptap(a_csr,desc_a,nlaggr,parms,ac,& nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) @@ -106,542 +104,38 @@ subroutine amg_d_ptap(a_csr,desc_a,nlaggr,parms,ac,& & desc_ac%get_local_rows(),desc_a%get_local_cols() if (debug) flush(0) - if (do_timings) call psb_tic(idx_spspmm) call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) if (debug) write(0,*) me,trim(name),' Done AxPROL ',& & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& & desc_ac%get_local_rows(),desc_ac%get_local_cols() - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - call coo_prol%cp_to_coo(coo_restr,info) - call coo_restr%set_ncols(desc_ac%get_local_cols()) - call coo_restr%set_nrows(desc_a%get_local_rows()) - call psb_d_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_d_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_d_ptap - -subroutine amg_d_ld_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_d_inner_mod - use amg_d_base_aggregator_mod !, amg_protect_name => amg_d_ld_ptap - implicit none - - ! Arguments - type(psb_d_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_dml_parms), intent(inout) :: parms - type(psb_ld_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_ldspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_d_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_ifmt(csr_prol,info) + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - type(psb_d_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_d_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_lcoo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_lcoo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_ld_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_d_ld_ptap - -subroutine amg_ld_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_d_inner_mod - use amg_d_base_aggregator_mod!, amg_protect_name => amg_ld_ptap - implicit none - - ! Arguments - type(psb_ld_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_dml_parms), intent(inout) :: parms - type(psb_ld_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_ldspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_ld_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr - - ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_fmt(csr_prol,info) - - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. + & 'starting sphalo/ rwxtd' - if (present(desc_ax)) then - block - type(psb_d_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_d_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) - write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& - & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() - if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) call ac_csr%set_nrows(desc_ac%get_local_rows()) call ac_csr%set_ncols(desc_ac%get_local_cols()) @@ -653,12 +147,8 @@ subroutine amg_ld_ptap(a_csr,desc_a,nlaggr,parms,ac,& ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() call coo_prol%set_ncols(desc_ac%get_local_cols()) - !call coo_restr%mv_from_ifmt(csr_restr,info) -!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) -!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done ptap ' @@ -669,21 +159,4 @@ subroutine amg_ld_ptap(a_csr,desc_a,nlaggr,parms,ac,& 9999 call psb_error_handler(err_act) return - -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_ld_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_ld_ptap +end subroutine amg_d_ptap diff --git a/mlprec/impl/aggregator/amg_d_ptap_bld.f90 b/mlprec/impl/aggregator/amg_d_ptap_bld.f90 new file mode 100644 index 00000000..0a29ca6b --- /dev/null +++ b/mlprec/impl/aggregator/amg_d_ptap_bld.f90 @@ -0,0 +1,697 @@ +! +! +! MLD2P4 Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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: amg_ptap_bld.F90 +! +! This routine does two things: +! 1. Computes the transpose of the prolongator Prol and stores +! it explicitly into Restr +! 2. Computes AC = P^T A P +! +! Yes, doing two things at once is against received wisdom, but we are +! trying to save time in the build procedure; maybe we will revisit this +! decision. +! +subroutine amg_d_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_d_inner_mod + use amg_d_base_aggregator_mod, amg_protect_name => amg_d_ptap_bld + implicit none + + ! Arguments + type(psb_d_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_d_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_dspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_d_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 + integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + call coo_prol%cp_to_coo(coo_restr,info) + call coo_restr%set_ncols(desc_ac%get_local_cols()) + call coo_restr%set_nrows(desc_a%get_local_rows()) + call psb_d_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_d_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_d_ptap_bld + +subroutine amg_d_ld_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_d_inner_mod + use amg_d_base_aggregator_mod !, amg_protect_name => amg_d_ld_ptap_bld + implicit none + + ! Arguments + type(psb_d_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_ld_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_ldspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_d_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_ifmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_d_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_d_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_lcoo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_lcoo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_ld_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_d_ld_ptap_bld + +subroutine amg_ld_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_d_inner_mod + use amg_d_base_aggregator_mod!, amg_protect_name => amg_ld_ptap_bld + implicit none + + ! Arguments + type(psb_ld_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_ld_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_ldspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_ld_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_d_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_d_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) + write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& + & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() + if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + !call coo_restr%mv_from_ifmt(csr_restr,info) +!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) +!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_ld_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_ld_ptap_bld diff --git a/mlprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 b/mlprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 index 3aba067c..8be8a25d 100644 --- a/mlprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 +++ b/mlprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 @@ -160,7 +160,7 @@ subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& call psb_cdasb(desc_ac,info) call psb_cd_reinit(desc_ac,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call coo_restr%set_nrows(desc_ac%get_local_rows()) diff --git a/mlprec/impl/aggregator/amg_daggrmat_smth_bld.f90 b/mlprec/impl/aggregator/amg_daggrmat_smth_bld.f90 index 1eb6ee83..e486cc68 100644 --- a/mlprec/impl/aggregator/amg_daggrmat_smth_bld.f90 +++ b/mlprec/impl/aggregator/amg_daggrmat_smth_bld.f90 @@ -286,7 +286,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& nzl = acsr1%get_nzeros() call acsr1%mv_to_coo(coo_prol,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call op_prol%mv_from(coo_prol) diff --git a/mlprec/impl/aggregator/amg_s_ptap.f90 b/mlprec/impl/aggregator/amg_s_ptap.f90 index 0270541d..ce6d71f1 100644 --- a/mlprec/impl/aggregator/amg_s_ptap.f90 +++ b/mlprec/impl/aggregator/amg_s_ptap.f90 @@ -32,11 +32,14 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -! File: amg_daggrmat_nosmth_bld.F90 +! File: amg_ptap.F90 ! +! This routine computes the triple product : +! AC = P^T A P +! P^T is stored in Restr ! subroutine amg_s_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) + & coo_prol,desc_ac,coo_restr,info) use psb_base_mod use amg_s_inner_mod use amg_s_base_aggregator_mod, amg_protect_name => amg_s_ptap @@ -51,7 +54,6 @@ subroutine amg_s_ptap(a_csr,desc_a,nlaggr,parms,ac,& type(psb_desc_type), intent(inout) :: desc_ac type(psb_sspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax ! Local variables integer(psb_ipk_) :: err_act @@ -64,8 +66,7 @@ subroutine amg_s_ptap(a_csr,desc_a,nlaggr,parms,ac,& integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 + logical, parameter :: debug=.false. name='amg_ptap' if(psb_get_errstatus().ne.0) return @@ -82,9 +83,6 @@ subroutine amg_s_ptap(a_csr,desc_a,nlaggr,parms,ac,& nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) @@ -106,542 +104,38 @@ subroutine amg_s_ptap(a_csr,desc_a,nlaggr,parms,ac,& & desc_ac%get_local_rows(),desc_a%get_local_cols() if (debug) flush(0) - if (do_timings) call psb_tic(idx_spspmm) call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) if (debug) write(0,*) me,trim(name),' Done AxPROL ',& & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& & desc_ac%get_local_rows(),desc_ac%get_local_cols() - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - call coo_prol%cp_to_coo(coo_restr,info) - call coo_restr%set_ncols(desc_ac%get_local_cols()) - call coo_restr%set_nrows(desc_a%get_local_rows()) - call psb_s_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_s_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_s_ptap - -subroutine amg_s_ls_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_s_inner_mod - use amg_s_base_aggregator_mod !, amg_protect_name => amg_s_ls_ptap - implicit none - - ! Arguments - type(psb_s_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_sml_parms), intent(inout) :: parms - type(psb_ls_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_lsspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_s_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_ifmt(csr_prol,info) + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - type(psb_s_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_s_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_lcoo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_lcoo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_ls_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_s_ls_ptap - -subroutine amg_ls_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_s_inner_mod - use amg_s_base_aggregator_mod!, amg_protect_name => amg_ls_ptap - implicit none - - ! Arguments - type(psb_ls_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_sml_parms), intent(inout) :: parms - type(psb_ls_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_lsspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_ls_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr - - ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_fmt(csr_prol,info) - - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. + & 'starting sphalo/ rwxtd' - if (present(desc_ax)) then - block - type(psb_s_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_s_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) - write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& - & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() - if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) call ac_csr%set_nrows(desc_ac%get_local_rows()) call ac_csr%set_ncols(desc_ac%get_local_cols()) @@ -653,12 +147,8 @@ subroutine amg_ls_ptap(a_csr,desc_a,nlaggr,parms,ac,& ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() call coo_prol%set_ncols(desc_ac%get_local_cols()) - !call coo_restr%mv_from_ifmt(csr_restr,info) -!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) -!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done ptap ' @@ -669,21 +159,4 @@ subroutine amg_ls_ptap(a_csr,desc_a,nlaggr,parms,ac,& 9999 call psb_error_handler(err_act) return - -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_ls_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_ls_ptap +end subroutine amg_s_ptap diff --git a/mlprec/impl/aggregator/amg_s_ptap_bld.f90 b/mlprec/impl/aggregator/amg_s_ptap_bld.f90 new file mode 100644 index 00000000..a3dad635 --- /dev/null +++ b/mlprec/impl/aggregator/amg_s_ptap_bld.f90 @@ -0,0 +1,697 @@ +! +! +! MLD2P4 Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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: amg_ptap_bld.F90 +! +! This routine does two things: +! 1. Computes the transpose of the prolongator Prol and stores +! it explicitly into Restr +! 2. Computes AC = P^T A P +! +! Yes, doing two things at once is against received wisdom, but we are +! trying to save time in the build procedure; maybe we will revisit this +! decision. +! +subroutine amg_s_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_s_inner_mod + use amg_s_base_aggregator_mod, amg_protect_name => amg_s_ptap_bld + implicit none + + ! Arguments + type(psb_s_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_sspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_s_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 + integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + call coo_prol%cp_to_coo(coo_restr,info) + call coo_restr%set_ncols(desc_ac%get_local_cols()) + call coo_restr%set_nrows(desc_a%get_local_rows()) + call psb_s_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_s_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_s_ptap_bld + +subroutine amg_s_ls_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_s_inner_mod + use amg_s_base_aggregator_mod !, amg_protect_name => amg_s_ls_ptap_bld + implicit none + + ! Arguments + type(psb_s_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_ls_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_lsspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_s_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_ifmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_s_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_s_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_lcoo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_lcoo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_ls_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_s_ls_ptap_bld + +subroutine amg_ls_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_s_inner_mod + use amg_s_base_aggregator_mod!, amg_protect_name => amg_ls_ptap_bld + implicit none + + ! Arguments + type(psb_ls_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_sml_parms), intent(inout) :: parms + type(psb_ls_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_lsspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_ls_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_s_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_s_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) + write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& + & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() + if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + !call coo_restr%mv_from_ifmt(csr_restr,info) +!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) +!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_ls_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_ls_ptap_bld diff --git a/mlprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 b/mlprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 index 5b5321d2..068b1b4c 100644 --- a/mlprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 +++ b/mlprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 @@ -160,7 +160,7 @@ subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& call psb_cdasb(desc_ac,info) call psb_cd_reinit(desc_ac,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call coo_restr%set_nrows(desc_ac%get_local_rows()) diff --git a/mlprec/impl/aggregator/amg_saggrmat_smth_bld.f90 b/mlprec/impl/aggregator/amg_saggrmat_smth_bld.f90 index eba0eda6..e74ba228 100644 --- a/mlprec/impl/aggregator/amg_saggrmat_smth_bld.f90 +++ b/mlprec/impl/aggregator/amg_saggrmat_smth_bld.f90 @@ -286,7 +286,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& nzl = acsr1%get_nzeros() call acsr1%mv_to_coo(coo_prol,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call op_prol%mv_from(coo_prol) diff --git a/mlprec/impl/aggregator/amg_z_ptap.f90 b/mlprec/impl/aggregator/amg_z_ptap.f90 index bea71826..6d299f69 100644 --- a/mlprec/impl/aggregator/amg_z_ptap.f90 +++ b/mlprec/impl/aggregator/amg_z_ptap.f90 @@ -32,11 +32,14 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! -! File: amg_daggrmat_nosmth_bld.F90 +! File: amg_ptap.F90 ! +! This routine computes the triple product : +! AC = P^T A P +! P^T is stored in Restr ! subroutine amg_z_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) + & coo_prol,desc_ac,coo_restr,info) use psb_base_mod use amg_z_inner_mod use amg_z_base_aggregator_mod, amg_protect_name => amg_z_ptap @@ -51,7 +54,6 @@ subroutine amg_z_ptap(a_csr,desc_a,nlaggr,parms,ac,& type(psb_desc_type), intent(inout) :: desc_ac type(psb_zspmat_type), intent(out) :: ac integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax ! Local variables integer(psb_ipk_) :: err_act @@ -64,8 +66,7 @@ subroutine amg_z_ptap(a_csr,desc_a,nlaggr,parms,ac,& integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 + logical, parameter :: debug=.false. name='amg_ptap' if(psb_get_errstatus().ne.0) return @@ -82,9 +83,6 @@ subroutine amg_z_ptap(a_csr,desc_a,nlaggr,parms,ac,& nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) @@ -106,542 +104,38 @@ subroutine amg_z_ptap(a_csr,desc_a,nlaggr,parms,ac,& & desc_ac%get_local_rows(),desc_a%get_local_cols() if (debug) flush(0) - if (do_timings) call psb_tic(idx_spspmm) call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) if (debug) write(0,*) me,trim(name),' Done AxPROL ',& & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& & desc_ac%get_local_rows(),desc_ac%get_local_cols() - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - call coo_prol%cp_to_coo(coo_restr,info) - call coo_restr%set_ncols(desc_ac%get_local_cols()) - call coo_restr%set_nrows(desc_a%get_local_rows()) - call psb_z_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_z_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_z_ptap - -subroutine amg_z_lz_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_z_inner_mod - use amg_z_base_aggregator_mod !, amg_protect_name => amg_z_lz_ptap - implicit none - - ! Arguments - type(psb_z_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_dml_parms), intent(inout) :: parms - type(psb_lz_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_lzspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_z_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_ifmt(csr_prol,info) + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. - - if (present(desc_ax)) then - block - type(psb_z_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_z_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) -! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& -! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() -! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_lcoo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_lcoo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) - - call ac_csr%set_nrows(desc_ac%get_local_rows()) - call ac_csr%set_ncols(desc_ac%get_local_cols()) - call ac%mv_from(ac_csr) - call ac%set_asb() - - if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() - if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr - ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() - - call coo_prol%set_ncols(desc_ac%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done ptap ' - - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(err_act) - - return -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_lz_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_z_lz_ptap - -subroutine amg_lz_ptap(a_csr,desc_a,nlaggr,parms,ac,& - & coo_prol,desc_ac,coo_restr,info,desc_ax) - use psb_base_mod - use amg_z_inner_mod - use amg_z_base_aggregator_mod!, amg_protect_name => amg_lz_ptap - implicit none - - ! Arguments - type(psb_lz_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(inout) :: desc_a - integer(psb_lpk_), intent(inout) :: nlaggr(:) - type(amg_dml_parms), intent(inout) :: parms - type(psb_lz_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr - type(psb_desc_type), intent(inout) :: desc_ac - type(psb_lzspmat_type), intent(out) :: ac - integer(psb_ipk_), intent(out) :: info - type(psb_desc_type), intent(inout), optional :: desc_ax - - ! Local variables - integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo - character(len=40) :: name - integer(psb_ipk_) :: ierr(5) - type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo - type(psb_lz_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr - integer(psb_ipk_) :: debug_level, debug_unit, naggr - integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & - & nzt, naggrm1, naggrp1, i, k - integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza - logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. - integer(psb_ipk_), save :: idx_spspmm=-1 - - name='amg_ptap' - if(psb_get_errstatus().ne.0) return - info=psb_success_ - call psb_erractionsave(err_act) - - - ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() - call psb_info(ictxt, me, np) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - nglob = desc_a%get_global_rows() - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() - - if ((do_timings).and.(idx_spspmm==-1)) & - & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") - - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr - - ! - ! COO_PROL should arrive here with local numbering - ! - if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& - & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& - & nrow,ntaggr,naggr - - call coo_prol%cp_to_fmt(csr_prol,info) - - if (debug) write(0,*) me,trim(name),' Product AxPROL ',& - & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & - & desc_a%get_local_rows(),desc_a%get_local_cols(),& - & desc_ac%get_local_rows(),desc_a%get_local_cols() - if (debug) flush(0) - - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - - if (debug) write(0,*) me,trim(name),' Done AxPROL ',& - & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& - & desc_ac%get_local_rows(),desc_ac%get_local_cols() - - ! - ! Ok first product done. + & 'starting sphalo/ rwxtd' - if (present(desc_ax)) then - block - type(psb_z_coo_sparse_mat) :: icoo_restr - - call coo_prol%cp_to_icoo(icoo_restr,info) - call icoo_restr%set_ncols(desc_ac%get_local_cols()) - call icoo_restr%set_nrows(desc_a%get_local_rows()) - call psb_z_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) - call icoo_restr%set_nrows(desc_ac%get_local_rows()) - call icoo_restr%set_ncols(desc_ax%get_local_cols()) - write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& - & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() - if (desc_a%get_local_cols()= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - - else - - ! - ! Remember that RESTR must be built from PROL after halo extension, - ! which is done above in psb_par_spspmm - if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& - & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() - call csr_prol%mv_to_coo(coo_restr,info) -!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& -!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() - if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) - - call coo_restr%transp() - nzl = coo_restr%get_nzeros() - nrl = desc_ac%get_local_rows() - i=0 - ! - ! Only keep local rows - ! - do k=1, nzl - if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then - i = i+1 - coo_restr%val(i) = coo_restr%val(k) - coo_restr%ia(i) = coo_restr%ia(k) - coo_restr%ja(i) = coo_restr%ja(k) - end if - end do - call coo_restr%set_nzeros(i) - call coo_restr%fix(info) - nzl = coo_restr%get_nzeros() - call coo_restr%set_nrows(desc_ac%get_local_rows()) - call coo_restr%set_ncols(desc_a%get_local_cols()) - if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) - call csr_restr%cp_from_coo(coo_restr,info) - -!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& - & csr_restr%get_nrows(),csr_restr%get_ncols(), & - & desc_ac%get_local_rows(),desc_a%get_local_cols(),& - & acsr3%get_nrows(),acsr3%get_ncols() - if (do_timings) call psb_tic(idx_spspmm) - call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) - if (do_timings) call psb_toc(idx_spspmm) - call acsr3%free() - end if - - if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) call ac_csr%set_nrows(desc_ac%get_local_rows()) call ac_csr%set_ncols(desc_ac%get_local_cols()) @@ -653,12 +147,8 @@ subroutine amg_lz_ptap(a_csr,desc_a,nlaggr,parms,ac,& ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() call coo_prol%set_ncols(desc_ac%get_local_cols()) - !call coo_restr%mv_from_ifmt(csr_restr,info) -!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) -!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done ptap ' @@ -669,21 +159,4 @@ subroutine amg_lz_ptap(a_csr,desc_a,nlaggr,parms,ac,& 9999 call psb_error_handler(err_act) return - -contains - subroutine check_coo(me,string,coo) - implicit none - integer(psb_ipk_) :: me - type(psb_lz_coo_sparse_mat) :: coo - character(len=*) :: string - integer(psb_lpk_) :: nr,nc,nz - nr = coo%get_nrows() - nc = coo%get_ncols() - nz = coo%get_nzeros() - write(0,*) me,string,nr,nc,& - & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& - & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) - - end subroutine check_coo - -end subroutine amg_lz_ptap +end subroutine amg_z_ptap diff --git a/mlprec/impl/aggregator/amg_z_ptap_bld.f90 b/mlprec/impl/aggregator/amg_z_ptap_bld.f90 new file mode 100644 index 00000000..54696ecf --- /dev/null +++ b/mlprec/impl/aggregator/amg_z_ptap_bld.f90 @@ -0,0 +1,697 @@ +! +! +! MLD2P4 Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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: amg_ptap_bld.F90 +! +! This routine does two things: +! 1. Computes the transpose of the prolongator Prol and stores +! it explicitly into Restr +! 2. Computes AC = P^T A P +! +! Yes, doing two things at once is against received wisdom, but we are +! trying to save time in the build procedure; maybe we will revisit this +! decision. +! +subroutine amg_z_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_z_inner_mod + use amg_z_base_aggregator_mod, amg_protect_name => amg_z_ptap_bld + implicit none + + ! Arguments + type(psb_z_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_zspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_z_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nglob, ntaggr, naggrm1, naggrp1 + integer(psb_ipk_) :: nrow, ncol, nrl, nzl, ip, nzt, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + call coo_prol%cp_to_coo(coo_restr,info) + call coo_restr%set_ncols(desc_ac%get_local_cols()) + call coo_restr%set_nrows(desc_a%get_local_rows()) + call psb_z_coo_glob_transpose(coo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),coo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_z_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_z_ptap_bld + +subroutine amg_z_lz_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_z_inner_mod + use amg_z_base_aggregator_mod !, amg_protect_name => amg_z_lz_ptap_bld + implicit none + + ! Arguments + type(psb_z_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_lz_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_lzspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_z_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_ifmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_z_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_z_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) +! !$ write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& +! !$ & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() +! !$ if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_lcoo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_lcoo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_lz_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_z_lz_ptap_bld + +subroutine amg_lz_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info,desc_ax) + use psb_base_mod + use amg_z_inner_mod + use amg_z_base_aggregator_mod!, amg_protect_name => amg_lz_ptap_bld + implicit none + + ! Arguments + type(psb_lz_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_lz_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr + type(psb_desc_type), intent(inout) :: desc_ac + type(psb_lzspmat_type), intent(out) :: ac + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(inout), optional :: desc_ax + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo + character(len=40) :: name + integer(psb_ipk_) :: ierr(5) + type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo + type(psb_lz_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza + logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1 + + name='amg_ptap_bld' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm") + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + !write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr + + ! + ! COO_PROL should arrive here with local numbering + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call coo_prol%cp_to_fmt(csr_prol,info) + + if (debug) write(0,*) me,trim(name),' Product AxPROL ',& + & a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), & + & desc_a%get_local_rows(),desc_a%get_local_cols(),& + & desc_ac%get_local_rows(),desc_a%get_local_cols() + if (debug) flush(0) + + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + + if (debug) write(0,*) me,trim(name),' Done AxPROL ',& + & acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),& + & desc_ac%get_local_rows(),desc_ac%get_local_cols() + + ! + ! Ok first product done. + + if (present(desc_ax)) then + block + type(psb_z_coo_sparse_mat) :: icoo_restr + + call coo_prol%cp_to_icoo(icoo_restr,info) + call icoo_restr%set_ncols(desc_ac%get_local_cols()) + call icoo_restr%set_nrows(desc_a%get_local_rows()) + call psb_z_coo_glob_transpose(icoo_restr,desc_a,info,desc_c=desc_ac,desc_rx=desc_ax) + call icoo_restr%set_nrows(desc_ac%get_local_rows()) + call icoo_restr%set_ncols(desc_ax%get_local_cols()) + write(0,*) me,' ',trim(name),' check on glob_transpose 1: ',& + & desc_a%get_local_cols(),desc_ax%get_local_cols(),icoo_restr%get_nzeros() + if (desc_a%get_local_cols()= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_ax,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + + else + + ! + ! Remember that RESTR must be built from PROL after halo extension, + ! which is done above in psb_par_spspmm + if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',& + & csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros() + call csr_prol%mv_to_coo(coo_restr,info) +!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),& +!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros() + if (debug) call check_coo(me,trim(name)//' Check 1 (before transp) on coo_restr:',coo_restr) + + call coo_restr%transp() + nzl = coo_restr%get_nzeros() + nrl = desc_ac%get_local_rows() + i=0 + ! + ! Only keep local rows + ! + do k=1, nzl + if ((1 <= coo_restr%ia(k)) .and.(coo_restr%ia(k) <= nrl)) then + i = i+1 + coo_restr%val(i) = coo_restr%val(k) + coo_restr%ia(i) = coo_restr%ia(k) + coo_restr%ja(i) = coo_restr%ja(k) + end if + end do + call coo_restr%set_nzeros(i) + call coo_restr%fix(info) + nzl = coo_restr%get_nzeros() + call coo_restr%set_nrows(desc_ac%get_local_rows()) + call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 2 on coo_restr:',coo_restr) + call csr_restr%cp_from_coo(coo_restr,info) + +!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting sphalo/ rwxtd' + + if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',& + & csr_restr%get_nrows(),csr_restr%get_ncols(), & + & desc_ac%get_local_rows(),desc_a%get_local_cols(),& + & acsr3%get_nrows(),acsr3%get_ncols() + if (do_timings) call psb_tic(idx_spspmm) + call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_ac,info) + if (do_timings) call psb_toc(idx_spspmm) + call acsr3%free() + end if + + if (.not.desc_ac%is_asb()) call psb_cdasb(desc_ac,info) + + call ac_csr%set_nrows(desc_ac%get_local_rows()) + call ac_csr%set_ncols(desc_ac%get_local_cols()) + call ac%mv_from(ac_csr) + call ac%set_asb() + + if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus() + if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr + ! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros() + + call coo_prol%set_ncols(desc_ac%get_local_cols()) + !call coo_restr%mv_from_ifmt(csr_restr,info) +!!$ call coo_restr%set_nrows(desc_ac%get_local_rows()) +!!$ call coo_restr%set_ncols(desc_a%get_local_cols()) + if (debug) call check_coo(me,trim(name)//' Check 3 on coo_restr:',coo_restr) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ptap ' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + subroutine check_coo(me,string,coo) + implicit none + integer(psb_ipk_) :: me + type(psb_lz_coo_sparse_mat) :: coo + character(len=*) :: string + integer(psb_lpk_) :: nr,nc,nz + nr = coo%get_nrows() + nc = coo%get_ncols() + nz = coo%get_nzeros() + write(0,*) me,string,nr,nc,& + & minval(coo%ia(1:nz)),maxval(coo%ia(1:nz)),& + & minval(coo%ja(1:nz)),maxval(coo%ja(1:nz)) + + end subroutine check_coo + +end subroutine amg_lz_ptap_bld diff --git a/mlprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 b/mlprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 index 8ff0f7e5..271e1dd7 100644 --- a/mlprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 +++ b/mlprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 @@ -160,7 +160,7 @@ subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& call psb_cdasb(desc_ac,info) call psb_cd_reinit(desc_ac,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call coo_restr%set_nrows(desc_ac%get_local_rows()) diff --git a/mlprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 b/mlprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 index 558105ea..3cfc65fe 100644 --- a/mlprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 +++ b/mlprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 @@ -286,7 +286,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& nzl = acsr1%get_nzeros() call acsr1%mv_to_coo(coo_prol,info) - call amg_ptap(acsr,desc_a,nlaggr,parms,ac,& + call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,& & coo_prol,desc_ac,coo_restr,info) call op_prol%mv_from(coo_prol)