From 8966ecb4a6f7e9031c1a56a164a21ad810bbef03 Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Wed, 9 Oct 2024 09:43:43 +0000 Subject: [PATCH] First implementation of l1-aggregation --- amgprec/amg_base_prec_type.F90 | 15 ++++--- amgprec/amg_c_inner_mod.f90 | 3 +- amgprec/amg_d_inner_mod.f90 | 3 +- amgprec/amg_d_parmatch_aggregator_mod.F90 | 6 ++- amgprec/amg_s_inner_mod.f90 | 3 +- amgprec/amg_s_parmatch_aggregator_mod.F90 | 6 ++- amgprec/amg_z_inner_mod.f90 | 3 +- .../amg_c_dec_aggregator_mat_bld.f90 | 17 ++++---- .../aggregator/amg_caggrmat_minnrg_bld.f90 | 13 +++++- .../aggregator/amg_caggrmat_nosmth_bld.f90 | 14 +++++-- .../impl/aggregator/amg_caggrmat_smth_bld.f90 | 40 ++++++++++++++++--- .../amg_d_dec_aggregator_mat_bld.f90 | 17 ++++---- .../amg_d_parmatch_aggregator_mat_bld.F90 | 17 ++++---- .../aggregator/amg_d_parmatch_smth_bld.F90 | 33 +++++++++++++-- .../aggregator/amg_d_parmatch_unsmth_bld.F90 | 12 +++++- .../aggregator/amg_daggrmat_minnrg_bld.f90 | 13 +++++- .../aggregator/amg_daggrmat_nosmth_bld.f90 | 14 +++++-- .../impl/aggregator/amg_daggrmat_smth_bld.f90 | 40 ++++++++++++++++--- .../amg_s_dec_aggregator_mat_bld.f90 | 17 ++++---- .../amg_s_parmatch_aggregator_mat_bld.F90 | 17 ++++---- .../aggregator/amg_s_parmatch_smth_bld.F90 | 33 +++++++++++++-- .../aggregator/amg_s_parmatch_unsmth_bld.F90 | 12 +++++- .../aggregator/amg_saggrmat_minnrg_bld.f90 | 13 +++++- .../aggregator/amg_saggrmat_nosmth_bld.f90 | 14 +++++-- .../impl/aggregator/amg_saggrmat_smth_bld.f90 | 40 ++++++++++++++++--- .../amg_z_dec_aggregator_mat_bld.f90 | 17 ++++---- .../aggregator/amg_zaggrmat_minnrg_bld.f90 | 13 +++++- .../aggregator/amg_zaggrmat_nosmth_bld.f90 | 14 +++++-- .../impl/aggregator/amg_zaggrmat_smth_bld.f90 | 40 ++++++++++++++++--- 29 files changed, 383 insertions(+), 116 deletions(-) diff --git a/amgprec/amg_base_prec_type.F90 b/amgprec/amg_base_prec_type.F90 index 896bf13b..d99ca566 100644 --- a/amgprec/amg_base_prec_type.F90 +++ b/amgprec/amg_base_prec_type.F90 @@ -288,11 +288,12 @@ module amg_base_prec_type ! ! Legal values for entry: amg_aggr_prol_ ! - integer(psb_ipk_), parameter :: amg_no_smooth_ = 0 - integer(psb_ipk_), parameter :: amg_smooth_prol_ = 1 - integer(psb_ipk_), parameter :: amg_min_energy_ = 2 + integer(psb_ipk_), parameter :: amg_no_smooth_ = 0 + integer(psb_ipk_), parameter :: amg_smooth_prol_ = 1 + integer(psb_ipk_), parameter :: amg_l1_smooth_prol_ = 2 + integer(psb_ipk_), parameter :: amg_min_energy_ = 3 ! Disabling min_energy for the time being. - integer(psb_ipk_), parameter :: amg_max_aggr_prol_=amg_smooth_prol_ + integer(psb_ipk_), parameter :: amg_max_aggr_prol_= amg_l1_smooth_prol_ ! ! Legal values for entry: amg_aggr_filter_ ! @@ -376,8 +377,8 @@ module amg_base_prec_type character(len=19), parameter, private :: & & eigen_estimates(0:0)=(/'infinity norm '/) character(len=15), parameter, private :: & - & aggr_prols(0:3)=(/'unsmoothed ','smoothed ',& - & 'min energy ','bizr. smoothed'/) + & aggr_prols(0:4)=(/'unsmoothed ','smoothed ',& + & 'l1-smoothed ','min energy ','bizr. smoothed'/) character(len=15), parameter, private :: & & aggr_filters(0:1)=(/'no filtering ','filtering '/) character(len=15), parameter, private :: & @@ -548,6 +549,8 @@ contains val = amg_no_smooth_ case('SMOOTHED') val = amg_smooth_prol_ + case('L1-SMOOTHED','L1SMOOTHED') + val = amg_l1_smooth_prol_ case('MINENERGY') val = amg_min_energy_ case('NOPREC') diff --git a/amgprec/amg_c_inner_mod.f90 b/amgprec/amg_c_inner_mod.f90 index ac260134..c97b8c3f 100644 --- a/amgprec/amg_c_inner_mod.f90 +++ b/amgprec/amg_c_inner_mod.f90 @@ -109,11 +109,12 @@ module amg_c_inner_mod end interface amg_map_to_tprol abstract interface - subroutine amg_caggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_caggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lcspmat_type import :: amg_c_onelev_type, amg_sml_parms implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) diff --git a/amgprec/amg_d_inner_mod.f90 b/amgprec/amg_d_inner_mod.f90 index 8fa96609..83176af7 100644 --- a/amgprec/amg_d_inner_mod.f90 +++ b/amgprec/amg_d_inner_mod.f90 @@ -109,11 +109,12 @@ module amg_d_inner_mod end interface amg_map_to_tprol abstract interface - subroutine amg_daggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_daggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_ldspmat_type import :: amg_d_onelev_type, amg_dml_parms implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) diff --git a/amgprec/amg_d_parmatch_aggregator_mod.F90 b/amgprec/amg_d_parmatch_aggregator_mod.F90 index 525bb0c3..def1fba1 100644 --- a/amgprec/amg_d_parmatch_aggregator_mod.F90 +++ b/amgprec/amg_d_parmatch_aggregator_mod.F90 @@ -244,11 +244,12 @@ module amg_d_parmatch_aggregator_mod end interface interface - subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_d_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: amg_d_parmatch_aggregator_type, psb_desc_type, psb_dspmat_type,& & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a @@ -262,11 +263,12 @@ module amg_d_parmatch_aggregator_mod end interface interface - subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: amg_d_parmatch_aggregator_type, psb_desc_type, psb_dspmat_type,& & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a diff --git a/amgprec/amg_s_inner_mod.f90 b/amgprec/amg_s_inner_mod.f90 index 85ed1089..be883b21 100644 --- a/amgprec/amg_s_inner_mod.f90 +++ b/amgprec/amg_s_inner_mod.f90 @@ -109,11 +109,12 @@ module amg_s_inner_mod end interface amg_map_to_tprol abstract interface - subroutine amg_saggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_saggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lsspmat_type import :: amg_s_onelev_type, amg_sml_parms implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) diff --git a/amgprec/amg_s_parmatch_aggregator_mod.F90 b/amgprec/amg_s_parmatch_aggregator_mod.F90 index d58bd750..059c5b73 100644 --- a/amgprec/amg_s_parmatch_aggregator_mod.F90 +++ b/amgprec/amg_s_parmatch_aggregator_mod.F90 @@ -244,11 +244,12 @@ module amg_s_parmatch_aggregator_mod end interface interface - subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_s_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: amg_s_parmatch_aggregator_type, psb_desc_type, psb_sspmat_type,& & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a @@ -262,11 +263,12 @@ module amg_s_parmatch_aggregator_mod end interface interface - subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_s_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: amg_s_parmatch_aggregator_type, psb_desc_type, psb_sspmat_type,& & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a diff --git a/amgprec/amg_z_inner_mod.f90 b/amgprec/amg_z_inner_mod.f90 index bf997651..fb7139fd 100644 --- a/amgprec/amg_z_inner_mod.f90 +++ b/amgprec/amg_z_inner_mod.f90 @@ -109,11 +109,12 @@ module amg_z_inner_mod end interface amg_map_to_tprol abstract interface - subroutine amg_zaggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,& + subroutine amg_zaggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,& & ac,desc_ac,op_prol,op_restr,t_prol,info) import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_lzspmat_type import :: amg_z_onelev_type, amg_dml_parms implicit none + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) diff --git a/amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90 b/amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90 index 2c9317d1..186e7a1e 100644 --- a/amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90 +++ b/amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90 @@ -177,23 +177,24 @@ subroutine amg_c_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& select case (parms%aggr_prol) case (amg_no_smooth_) - call amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,& - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_caggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - case(amg_smooth_prol_) + case(amg_smooth_prol_,amg_l1_smooth_prol_) - call amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_caggrmat_smth_bld(parms%aggr_prol,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,& + op_restr,t_prol,info) !!$ case(amg_biz_prol_) !!$ !!$ call amg_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, & !!$ & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - + case(amg_min_energy_) - call amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_caggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) case default info = psb_err_internal_error_ diff --git a/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 index 8ac86dc6..3d223279 100644 --- a/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 @@ -69,6 +69,7 @@ ! ! ! Arguments: +! dol1smoothing - fictitious integer argument, it is not used inside ! a - type(psb_cspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -104,8 +105,8 @@ ! Error code. ! ! -subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_c_inner_mod, amg_protect_name => amg_caggrmat_minnrg_bld @@ -113,6 +114,7 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -171,6 +173,13 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + + !NEEDS TO BE REWORKED !! ! naggr: number of local aggregates diff --git a/amgprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 b/amgprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 index 87c79dc6..9699545e 100644 --- a/amgprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 +++ b/amgprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90 @@ -94,10 +94,11 @@ ! ! info - integer, output. ! Error code. +! dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the +! code ! -! -subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_caggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_c_inner_mod, amg_protect_name => amg_caggrmat_nosmth_bld @@ -105,6 +106,7 @@ subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -137,6 +139,12 @@ subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& ctxt = desc_a%get_context() call psb_info(ctxt, me, np) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() diff --git a/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 index 67fce476..4cd75284 100644 --- a/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 @@ -69,6 +69,8 @@ ! ! ! Arguments: +! dol1smooth - Integer taking the type of smoother that has to be used +! on the tentative prolongator ! a - type(psb_cspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -102,8 +104,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_c_inner_mod, amg_protect_name => amg_caggrmat_smth_bld @@ -112,6 +114,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -132,7 +135,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_c_coo_sparse_mat) :: coo_prol, coo_restr type(psb_c_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr complex(psb_spk_), allocatable :: adiag(:) - real(psb_spk_), allocatable :: arwsum(:) + real(psb_spk_), allocatable :: arwsum(:),l1rwsum(:) integer(psb_ipk_) :: ierr(5) logical :: filter_mat integer(psb_ipk_) :: debug_level, debug_unit, err_act @@ -141,6 +144,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& logical, parameter :: debug_new=.false. character(len=80) :: filename logical, parameter :: do_timings=.false. + logical :: do_l1correction=.false. integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1 integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1 @@ -173,6 +177,9 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if ((do_timings).and.(idx_ptap==-1)) & & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld ") + ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator + if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true. + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() @@ -200,6 +207,24 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (info == psb_success_) & & call psb_halo(adiag,desc_a,info) if (info == psb_success_) call a%cp_to(acsr) + ! + ! Do the l1-correction on the diagonal if it is requested + ! + if (do_l1correction) then + allocate(l1rwsum(nrow)) + call acsr%arwsum(l1rwsum) + if (info == psb_success_) & + & call psb_realloc(ncol,l1rwsum,info) + if (info == psb_success_) & + & call psb_halo(l1rwsum,desc_a,info) + ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}| + !$OMP parallel do private(i) schedule(static) + do i=1,size(adiag) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) + end do + !$OMP end parallel do + end if + if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') @@ -230,7 +255,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& enddo if (jd == -1) then - write(0,*) name,': Warning: there is no diagonal element', i + if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i else acsrf%val(jd)=acsrf%val(jd)-tmp end if @@ -252,14 +277,16 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - if (parms%aggr_eig == amg_max_norm_) then + if (do_l1correction) then + ! For l1-Jacobi this can be estimated with 1 + parms%aggr_omega_val = done + else if (parms%aggr_eig == amg_max_norm_) then allocate(arwsum(nrow)) call acsr%arwsum(arwsum) anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow))) call psb_amx(ctxt,anorm) omega = 4.d0/(3.d0*anorm) parms%aggr_omega_val = omega - else info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid amg_aggr_eig_') @@ -322,6 +349,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done smooth_aggregate ' + if (allocated(l1rwsum)) deallocate(l1rwsum) call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90 b/amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90 index 7b01a0b8..65bd08e8 100644 --- a/amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90 +++ b/amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90 @@ -177,23 +177,24 @@ subroutine amg_d_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& select case (parms%aggr_prol) case (amg_no_smooth_) - call amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,& - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_daggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - case(amg_smooth_prol_) + case(amg_smooth_prol_,amg_l1_smooth_prol_) - call amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_daggrmat_smth_bld(parms%aggr_prol,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,& + op_restr,t_prol,info) !!$ case(amg_biz_prol_) !!$ !!$ call amg_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, & !!$ & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - + case(amg_min_energy_) - call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) case default info = psb_err_internal_error_ diff --git a/amgprec/impl/aggregator/amg_d_parmatch_aggregator_mat_bld.F90 b/amgprec/impl/aggregator/amg_d_parmatch_aggregator_mat_bld.F90 index 9b1171e0..1766d4a4 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_aggregator_mat_bld.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_aggregator_mat_bld.F90 @@ -184,20 +184,23 @@ subroutine amg_d_parmatch_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& ! select case (parms%aggr_prol) case (amg_no_smooth_) - call amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_d_parmatch_unsmth_bld(parms%aggr_prol,ag,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,& + t_prol,info) - case(amg_smooth_prol_) - call amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) + case(amg_smooth_prol_,amg_l1_smooth_prol_) + call amg_d_parmatch_smth_bld(parms%aggr_prol,ag,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,& + t_prol,info) !!$ case(amg_biz_prol_) !!$ call amg_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, & !!$ & parms,ac,desc_ac,op_prol,op_restr,info) case(amg_min_energy_) - call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,& + t_prol,info) case default info = psb_err_internal_error_ diff --git a/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 b/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 index b20b2fa8..23cd1459 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 @@ -69,6 +69,8 @@ ! ! ! Arguments: +! dol1smoothing - Select between l1-Jacobi and Jacobi as smoother for the +! tentative prolongator ! a - type(psb_dspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -102,8 +104,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_d_inner_mod @@ -116,6 +118,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a @@ -137,7 +140,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& type(psb_d_coo_sparse_mat) :: coo_prol, coo_restr type(psb_d_csr_sparse_mat) :: acsrf, csr_prol, acsr, tcsr real(psb_dpk_), allocatable :: adiag(:) - real(psb_dpk_), allocatable :: arwsum(:) + real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:) logical :: filter_mat integer(psb_ipk_) :: debug_level, debug_unit, err_act integer(psb_ipk_), parameter :: ncmax=16 @@ -145,6 +148,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& logical, parameter :: debug_new=.false., dump_r=.false., dump_p=.false., debug=.false. character(len=80) :: filename logical, parameter :: do_timings=.false. + logical :: do_l1correction=.false. integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1, idx_phase3=-1 integer(psb_ipk_), save :: idx_cdasb=-1, idx_ptap=-1 @@ -166,6 +170,10 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& ncol = desc_a%get_local_cols() theta = parms%aggr_thresh + ! Check if we have to perform l1-Jacobi or Jacobi as smoother + if(dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true. + + !write(0,*) me,' ',trim(name),' Start ',idx_spspmm if ((do_timings).and.(idx_spspmm==-1)) & & idx_spspmm = psb_get_timer_idx("PMC_SMTH_BLD: par_spspmm") @@ -217,6 +225,19 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& if (info == psb_success_) & & call psb_halo(adiag,desc_a,info) if (info == psb_success_) call a%cp_to(acsr) + ! Get the l1-diagonal of D + if (do_l1correction) then + allocate(l1rwsum(nrow)) + call acsr%arwsum(l1rwsum) + if (info == psb_success_) & + & call psb_realloc(ncol,l1rwsum,info) + if (info == psb_success_) & + & call psb_halo(l1rwsum,desc_a,info) + ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}| + do i=1,size(adiag) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) + end do + end if if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') @@ -267,7 +288,10 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& if (parms%aggr_omega_alg == amg_eig_est_) then - if (parms%aggr_eig == amg_max_norm_) then + if (do_l1correction) then + ! For l1-Jacobi this can be estimated with 1 + parms%aggr_omega_val = done + else if (parms%aggr_eig == amg_max_norm_) then allocate(arwsum(nrow)) call acsr%arwsum(arwsum) anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow))) @@ -373,6 +397,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& end block end if + if (allocated(l1rwsum)) deallocate(l1rwsum) if (do_timings) call psb_toc(idx_phase2) if (debug_level >= psb_debug_outer_) & diff --git a/amgprec/impl/aggregator/amg_d_parmatch_unsmth_bld.F90 b/amgprec/impl/aggregator/amg_d_parmatch_unsmth_bld.F90 index dc6574a0..85f1ec28 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_unsmth_bld.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_unsmth_bld.F90 @@ -68,6 +68,8 @@ ! ! ! Arguments: +! dol1smoothing - this not actually used inside unsmoothed aggregation, it +! is used just to perform a check ! a - type(psb_dspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -101,8 +103,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_d_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_d_inner_mod @@ -115,6 +117,7 @@ subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a @@ -159,6 +162,11 @@ subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& ictxt = desc_a%get_context() call psb_info(ictxt, me, np) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if #if !defined(SERIAL_MPI) nglob = desc_a%get_global_rows() diff --git a/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 index 510339e7..c9f87f72 100644 --- a/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 @@ -69,6 +69,7 @@ ! ! ! Arguments: +! dol1smoothing - fictitious integer argument, it is not used inside ! a - type(psb_dspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -104,8 +105,8 @@ ! Error code. ! ! -subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_d_inner_mod, amg_protect_name => amg_daggrmat_minnrg_bld @@ -113,6 +114,7 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -171,6 +173,13 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + + !NEEDS TO BE REWORKED !! ! naggr: number of local aggregates diff --git a/amgprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 b/amgprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 index 78e396cc..345316b1 100644 --- a/amgprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 +++ b/amgprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90 @@ -94,10 +94,11 @@ ! ! info - integer, output. ! Error code. +! dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the +! code ! -! -subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_daggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_d_inner_mod, amg_protect_name => amg_daggrmat_nosmth_bld @@ -105,6 +106,7 @@ subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -137,6 +139,12 @@ subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& ctxt = desc_a%get_context() call psb_info(ctxt, me, np) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() diff --git a/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 index d12f9c6b..399e27b6 100644 --- a/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 @@ -69,6 +69,8 @@ ! ! ! Arguments: +! dol1smooth - Integer taking the type of smoother that has to be used +! on the tentative prolongator ! a - type(psb_dspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -102,8 +104,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_d_inner_mod, amg_protect_name => amg_daggrmat_smth_bld @@ -112,6 +114,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -132,7 +135,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_d_coo_sparse_mat) :: coo_prol, coo_restr type(psb_d_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr real(psb_dpk_), allocatable :: adiag(:) - real(psb_dpk_), allocatable :: arwsum(:) + real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:) integer(psb_ipk_) :: ierr(5) logical :: filter_mat integer(psb_ipk_) :: debug_level, debug_unit, err_act @@ -141,6 +144,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& logical, parameter :: debug_new=.false. character(len=80) :: filename logical, parameter :: do_timings=.false. + logical :: do_l1correction=.false. integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1 integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1 @@ -173,6 +177,9 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if ((do_timings).and.(idx_ptap==-1)) & & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld ") + ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator + if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true. + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() @@ -200,6 +207,24 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (info == psb_success_) & & call psb_halo(adiag,desc_a,info) if (info == psb_success_) call a%cp_to(acsr) + ! + ! Do the l1-correction on the diagonal if it is requested + ! + if (do_l1correction) then + allocate(l1rwsum(nrow)) + call acsr%arwsum(l1rwsum) + if (info == psb_success_) & + & call psb_realloc(ncol,l1rwsum,info) + if (info == psb_success_) & + & call psb_halo(l1rwsum,desc_a,info) + ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}| + !$OMP parallel do private(i) schedule(static) + do i=1,size(adiag) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) + end do + !$OMP end parallel do + end if + if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') @@ -230,7 +255,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& enddo if (jd == -1) then - write(0,*) name,': Warning: there is no diagonal element', i + if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i else acsrf%val(jd)=acsrf%val(jd)-tmp end if @@ -252,14 +277,16 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - if (parms%aggr_eig == amg_max_norm_) then + if (do_l1correction) then + ! For l1-Jacobi this can be estimated with 1 + parms%aggr_omega_val = done + else if (parms%aggr_eig == amg_max_norm_) then allocate(arwsum(nrow)) call acsr%arwsum(arwsum) anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow))) call psb_amx(ctxt,anorm) omega = 4.d0/(3.d0*anorm) parms%aggr_omega_val = omega - else info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid amg_aggr_eig_') @@ -322,6 +349,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done smooth_aggregate ' + if (allocated(l1rwsum)) deallocate(l1rwsum) call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/aggregator/amg_s_dec_aggregator_mat_bld.f90 b/amgprec/impl/aggregator/amg_s_dec_aggregator_mat_bld.f90 index 39e96ec7..750889de 100644 --- a/amgprec/impl/aggregator/amg_s_dec_aggregator_mat_bld.f90 +++ b/amgprec/impl/aggregator/amg_s_dec_aggregator_mat_bld.f90 @@ -177,23 +177,24 @@ subroutine amg_s_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& select case (parms%aggr_prol) case (amg_no_smooth_) - call amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,& - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_saggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - case(amg_smooth_prol_) + case(amg_smooth_prol_,amg_l1_smooth_prol_) - call amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_saggrmat_smth_bld(parms%aggr_prol,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,& + op_restr,t_prol,info) !!$ case(amg_biz_prol_) !!$ !!$ call amg_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, & !!$ & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - + case(amg_min_energy_) - call amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_saggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) case default info = psb_err_internal_error_ diff --git a/amgprec/impl/aggregator/amg_s_parmatch_aggregator_mat_bld.F90 b/amgprec/impl/aggregator/amg_s_parmatch_aggregator_mat_bld.F90 index 8c10d006..38dff13c 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_aggregator_mat_bld.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_aggregator_mat_bld.F90 @@ -184,20 +184,23 @@ subroutine amg_s_parmatch_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& ! select case (parms%aggr_prol) case (amg_no_smooth_) - call amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_s_parmatch_unsmth_bld(parms%aggr_prol,ag,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,& + t_prol,info) - case(amg_smooth_prol_) - call amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) + case(amg_smooth_prol_,amg_l1_smooth_prol_) + call amg_s_parmatch_smth_bld(parms%aggr_prol,ag,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,& + t_prol,info) !!$ case(amg_biz_prol_) !!$ call amg_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, & !!$ & parms,ac,desc_ac,op_prol,op_restr,info) case(amg_min_energy_) - call amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_saggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,& + t_prol,info) case default info = psb_err_internal_error_ diff --git a/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 b/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 index 323853d1..7a53055f 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 @@ -69,6 +69,8 @@ ! ! ! Arguments: +! dol1smoothing - Select between l1-Jacobi and Jacobi as smoother for the +! tentative prolongator ! a - type(psb_sspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -102,8 +104,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_s_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_s_inner_mod @@ -116,6 +118,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a @@ -137,7 +140,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& type(psb_s_coo_sparse_mat) :: coo_prol, coo_restr type(psb_s_csr_sparse_mat) :: acsrf, csr_prol, acsr, tcsr real(psb_spk_), allocatable :: adiag(:) - real(psb_spk_), allocatable :: arwsum(:) + real(psb_spk_), allocatable :: arwsum(:),l1rwsum(:) logical :: filter_mat integer(psb_ipk_) :: debug_level, debug_unit, err_act integer(psb_ipk_), parameter :: ncmax=16 @@ -145,6 +148,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& logical, parameter :: debug_new=.false., dump_r=.false., dump_p=.false., debug=.false. character(len=80) :: filename logical, parameter :: do_timings=.false. + logical :: do_l1correction=.false. integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1, idx_phase3=-1 integer(psb_ipk_), save :: idx_cdasb=-1, idx_ptap=-1 @@ -166,6 +170,10 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& ncol = desc_a%get_local_cols() theta = parms%aggr_thresh + ! Check if we have to perform l1-Jacobi or Jacobi as smoother + if(dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true. + + !write(0,*) me,' ',trim(name),' Start ',idx_spspmm if ((do_timings).and.(idx_spspmm==-1)) & & idx_spspmm = psb_get_timer_idx("PMC_SMTH_BLD: par_spspmm") @@ -217,6 +225,19 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& if (info == psb_success_) & & call psb_halo(adiag,desc_a,info) if (info == psb_success_) call a%cp_to(acsr) + ! Get the l1-diagonal of D + if (do_l1correction) then + allocate(l1rwsum(nrow)) + call acsr%arwsum(l1rwsum) + if (info == psb_success_) & + & call psb_realloc(ncol,l1rwsum,info) + if (info == psb_success_) & + & call psb_halo(l1rwsum,desc_a,info) + ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}| + do i=1,size(adiag) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) + end do + end if if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') @@ -267,7 +288,10 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& if (parms%aggr_omega_alg == amg_eig_est_) then - if (parms%aggr_eig == amg_max_norm_) then + if (do_l1correction) then + ! For l1-Jacobi this can be estimated with 1 + parms%aggr_omega_val = done + else if (parms%aggr_eig == amg_max_norm_) then allocate(arwsum(nrow)) call acsr%arwsum(arwsum) anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow))) @@ -373,6 +397,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& end block end if + if (allocated(l1rwsum)) deallocate(l1rwsum) if (do_timings) call psb_toc(idx_phase2) if (debug_level >= psb_debug_outer_) & diff --git a/amgprec/impl/aggregator/amg_s_parmatch_unsmth_bld.F90 b/amgprec/impl/aggregator/amg_s_parmatch_unsmth_bld.F90 index 6fcf65ac..4a96ab95 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_unsmth_bld.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_unsmth_bld.F90 @@ -68,6 +68,8 @@ ! ! ! Arguments: +! dol1smoothing - this not actually used inside unsmoothed aggregation, it +! is used just to perform a check ! a - type(psb_sspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -101,8 +103,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_s_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_s_inner_mod @@ -115,6 +117,7 @@ subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a @@ -159,6 +162,11 @@ subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& ictxt = desc_a%get_context() call psb_info(ictxt, me, np) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if #if !defined(SERIAL_MPI) nglob = desc_a%get_global_rows() diff --git a/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 index a609e382..80a35344 100644 --- a/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 @@ -69,6 +69,7 @@ ! ! ! Arguments: +! dol1smoothing - fictitious integer argument, it is not used inside ! a - type(psb_sspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -104,8 +105,8 @@ ! Error code. ! ! -subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_saggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_s_inner_mod, amg_protect_name => amg_saggrmat_minnrg_bld @@ -113,6 +114,7 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -171,6 +173,13 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + + !NEEDS TO BE REWORKED !! ! naggr: number of local aggregates diff --git a/amgprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 b/amgprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 index ceeca998..fe684d0c 100644 --- a/amgprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 +++ b/amgprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90 @@ -94,10 +94,11 @@ ! ! info - integer, output. ! Error code. +! dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the +! code ! -! -subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_saggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_s_inner_mod, amg_protect_name => amg_saggrmat_nosmth_bld @@ -105,6 +106,7 @@ subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -137,6 +139,12 @@ subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& ctxt = desc_a%get_context() call psb_info(ctxt, me, np) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() diff --git a/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 index b70490ba..7cf657eb 100644 --- a/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 @@ -69,6 +69,8 @@ ! ! ! Arguments: +! dol1smooth - Integer taking the type of smoother that has to be used +! on the tentative prolongator ! a - type(psb_sspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -102,8 +104,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_s_inner_mod, amg_protect_name => amg_saggrmat_smth_bld @@ -112,6 +114,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -132,7 +135,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_s_coo_sparse_mat) :: coo_prol, coo_restr type(psb_s_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr real(psb_spk_), allocatable :: adiag(:) - real(psb_spk_), allocatable :: arwsum(:) + real(psb_spk_), allocatable :: arwsum(:),l1rwsum(:) integer(psb_ipk_) :: ierr(5) logical :: filter_mat integer(psb_ipk_) :: debug_level, debug_unit, err_act @@ -141,6 +144,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& logical, parameter :: debug_new=.false. character(len=80) :: filename logical, parameter :: do_timings=.false. + logical :: do_l1correction=.false. integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1 integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1 @@ -173,6 +177,9 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if ((do_timings).and.(idx_ptap==-1)) & & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld ") + ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator + if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true. + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() @@ -200,6 +207,24 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (info == psb_success_) & & call psb_halo(adiag,desc_a,info) if (info == psb_success_) call a%cp_to(acsr) + ! + ! Do the l1-correction on the diagonal if it is requested + ! + if (do_l1correction) then + allocate(l1rwsum(nrow)) + call acsr%arwsum(l1rwsum) + if (info == psb_success_) & + & call psb_realloc(ncol,l1rwsum,info) + if (info == psb_success_) & + & call psb_halo(l1rwsum,desc_a,info) + ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}| + !$OMP parallel do private(i) schedule(static) + do i=1,size(adiag) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) + end do + !$OMP end parallel do + end if + if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') @@ -230,7 +255,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& enddo if (jd == -1) then - write(0,*) name,': Warning: there is no diagonal element', i + if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i else acsrf%val(jd)=acsrf%val(jd)-tmp end if @@ -252,14 +277,16 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - if (parms%aggr_eig == amg_max_norm_) then + if (do_l1correction) then + ! For l1-Jacobi this can be estimated with 1 + parms%aggr_omega_val = done + else if (parms%aggr_eig == amg_max_norm_) then allocate(arwsum(nrow)) call acsr%arwsum(arwsum) anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow))) call psb_amx(ctxt,anorm) omega = 4.d0/(3.d0*anorm) parms%aggr_omega_val = omega - else info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid amg_aggr_eig_') @@ -322,6 +349,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done smooth_aggregate ' + if (allocated(l1rwsum)) deallocate(l1rwsum) call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/aggregator/amg_z_dec_aggregator_mat_bld.f90 b/amgprec/impl/aggregator/amg_z_dec_aggregator_mat_bld.f90 index 7135bfc3..e3b9f6af 100644 --- a/amgprec/impl/aggregator/amg_z_dec_aggregator_mat_bld.f90 +++ b/amgprec/impl/aggregator/amg_z_dec_aggregator_mat_bld.f90 @@ -177,23 +177,24 @@ subroutine amg_z_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& select case (parms%aggr_prol) case (amg_no_smooth_) - call amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,& - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_zaggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - case(amg_smooth_prol_) + case(amg_smooth_prol_,amg_l1_smooth_prol_) - call amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_zaggrmat_smth_bld(parms%aggr_prol,a,desc_a,& + ilaggr,nlaggr,parms,ac,desc_ac,op_prol,& + op_restr,t_prol,info) !!$ case(amg_biz_prol_) !!$ !!$ call amg_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, & !!$ & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) - + case(amg_min_energy_) - call amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, & - & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call amg_zaggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,& + nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info) case default info = psb_err_internal_error_ diff --git a/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 index 7bafbc18..eaa7f273 100644 --- a/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 @@ -69,6 +69,7 @@ ! ! ! Arguments: +! dol1smoothing - fictitious integer argument, it is not used inside ! a - type(psb_zspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -104,8 +105,8 @@ ! Error code. ! ! -subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_zaggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_minnrg_bld @@ -113,6 +114,7 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -171,6 +173,13 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + + !NEEDS TO BE REWORKED !! ! naggr: number of local aggregates diff --git a/amgprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 b/amgprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 index 6fd53861..2a0b631c 100644 --- a/amgprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 +++ b/amgprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90 @@ -94,10 +94,11 @@ ! ! info - integer, output. ! Error code. +! dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the +! code ! -! -subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_zaggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_nosmth_bld @@ -105,6 +106,7 @@ subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -137,6 +139,12 @@ subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,& ctxt = desc_a%get_context() call psb_info(ctxt, me, np) + if (dol1smoothing.ne.amg_no_smooth_) then + info=psb_err_fatal_; + call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?') + goto 9999 + end if + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() diff --git a/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 index 11f30589..70dd0ecd 100644 --- a/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 @@ -69,6 +69,8 @@ ! ! ! Arguments: +! dol1smooth - Integer taking the type of smoother that has to be used +! on the tentative prolongator ! a - type(psb_zspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. @@ -102,8 +104,8 @@ ! info - integer, output. ! Error code. ! -subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,t_prol,info) +subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& + parms,ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_smth_bld @@ -112,6 +114,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& implicit none ! Arguments + integer(psb_ipk_), intent(in) :: dol1smoothing type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) @@ -132,7 +135,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_z_coo_sparse_mat) :: coo_prol, coo_restr type(psb_z_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr complex(psb_dpk_), allocatable :: adiag(:) - real(psb_dpk_), allocatable :: arwsum(:) + real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:) integer(psb_ipk_) :: ierr(5) logical :: filter_mat integer(psb_ipk_) :: debug_level, debug_unit, err_act @@ -141,6 +144,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& logical, parameter :: debug_new=.false. character(len=80) :: filename logical, parameter :: do_timings=.false. + logical :: do_l1correction=.false. integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1 integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1 @@ -173,6 +177,9 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if ((do_timings).and.(idx_ptap==-1)) & & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld ") + ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator + if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true. + nglob = desc_a%get_global_rows() nrow = desc_a%get_local_rows() @@ -200,6 +207,24 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (info == psb_success_) & & call psb_halo(adiag,desc_a,info) if (info == psb_success_) call a%cp_to(acsr) + ! + ! Do the l1-correction on the diagonal if it is requested + ! + if (do_l1correction) then + allocate(l1rwsum(nrow)) + call acsr%arwsum(l1rwsum) + if (info == psb_success_) & + & call psb_realloc(ncol,l1rwsum,info) + if (info == psb_success_) & + & call psb_halo(l1rwsum,desc_a,info) + ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}| + !$OMP parallel do private(i) schedule(static) + do i=1,size(adiag) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) + end do + !$OMP end parallel do + end if + if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') @@ -230,7 +255,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& enddo if (jd == -1) then - write(0,*) name,': Warning: there is no diagonal element', i + if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i else acsrf%val(jd)=acsrf%val(jd)-tmp end if @@ -252,14 +277,16 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - if (parms%aggr_eig == amg_max_norm_) then + if (do_l1correction) then + ! For l1-Jacobi this can be estimated with 1 + parms%aggr_omega_val = done + else if (parms%aggr_eig == amg_max_norm_) then allocate(arwsum(nrow)) call acsr%arwsum(arwsum) anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow))) call psb_amx(ctxt,anorm) omega = 4.d0/(3.d0*anorm) parms%aggr_omega_val = omega - else info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid amg_aggr_eig_') @@ -322,6 +349,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done smooth_aggregate ' + if (allocated(l1rwsum)) deallocate(l1rwsum) call psb_erractionrestore(err_act) return