From 53afa87814afd74b4ed6721387f96f051963e43f Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Fri, 13 Dec 2024 17:19:04 +0000 Subject: [PATCH] Added row-preserving filtering --- amgprec/amg_base_prec_type.F90 | 12 ++++++---- .../impl/aggregator/amg_caggrmat_smth_bld.f90 | 23 +++++++++++-------- .../impl/aggregator/amg_daggrmat_smth_bld.f90 | 23 +++++++++++-------- .../impl/aggregator/amg_saggrmat_smth_bld.f90 | 23 +++++++++++-------- .../impl/aggregator/amg_zaggrmat_smth_bld.f90 | 23 +++++++++++-------- 5 files changed, 64 insertions(+), 40 deletions(-) diff --git a/amgprec/amg_base_prec_type.F90 b/amgprec/amg_base_prec_type.F90 index d99ca566..70e55801 100644 --- a/amgprec/amg_base_prec_type.F90 +++ b/amgprec/amg_base_prec_type.F90 @@ -297,9 +297,10 @@ module amg_base_prec_type ! ! Legal values for entry: amg_aggr_filter_ ! - integer(psb_ipk_), parameter :: amg_no_filter_mat_ = 0 - integer(psb_ipk_), parameter :: amg_filter_mat_ = 1 - integer(psb_ipk_), parameter :: amg_max_filter_mat_ = amg_filter_mat_ + integer(psb_ipk_), parameter :: amg_no_filter_mat_ = 0 + integer(psb_ipk_), parameter :: amg_filter_mat_ = 1 + integer(psb_ipk_), parameter :: amg_filter_prow_mat_ = 2 + integer(psb_ipk_), parameter :: amg_max_filter_mat_ = amg_filter_prow_mat_ ! ! Legal values for entry: amg_aggr_ord_ ! @@ -380,7 +381,8 @@ module amg_base_prec_type & aggr_prols(0:4)=(/'unsmoothed ','smoothed ',& & 'l1-smoothed ','min energy ','bizr. smoothed'/) character(len=15), parameter, private :: & - & aggr_filters(0:1)=(/'no filtering ','filtering '/) + & aggr_filters(0:2)=(/'no filtering ','filtering ',& + & 'filtering rsum'/) character(len=15), parameter, private :: & & matrix_names(0:1)=(/'distributed ','replicated '/) character(len=18), parameter, private :: & @@ -591,6 +593,8 @@ contains val = amg_eig_est_ case('FILTER') val = amg_filter_mat_ + case('FILTERROWSUM') + val = amg_filter_prow_mat_ case('NOFILTER','NO_FILTER') val = amg_no_filter_mat_ case('OUTER_SWEEPS') diff --git a/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 index e87637cf..209e4570 100644 --- a/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 @@ -110,6 +110,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& use amg_base_prec_type use amg_c_inner_mod, amg_protect_name => amg_caggrmat_smth_bld use amg_c_base_aggregator_mod +! use, intrinsic :: ieee_arithmetic implicit none @@ -192,7 +193,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (parms%aggr_filter == amg_filter_mat_) + filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_) ! ! naggr: number of local aggregates @@ -220,7 +221,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& ! \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)) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) end do !$OMP end parallel do end if @@ -257,8 +258,13 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& if (jd == -1) then ! if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i - else + else if (parms%aggr_filter == amg_filter_mat_) then + ! We perform filtering in the standard way assuming that A is an M-matrix acsrf%val(jd)=acsrf%val(jd)-tmp + else if (parms%aggr_filter == amg_filter_prow_mat_) then + ! We are probably doing l1-correction, hence we want to preserve the + ! row sum of the matrix: note the change in sign + acsrf%val(jd)=acsrf%val(jd)+tmp end if enddo !$OMP end parallel do @@ -266,7 +272,6 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& call acsrf%clean_zeros(info) end if - !$OMP parallel do private(i) schedule(static) do i=1,size(adiag) if (adiag(i) /= czero) then @@ -278,11 +283,11 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - !if (do_l1correction) then - ! ! For l1-Jacobi this can be estimated with 1 - ! parms%aggr_omega_val = done - ! - if (parms%aggr_eig == amg_max_norm_) then + if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then + ! For l1-Jacobi this can be estimated with 1: + ! this makes sense only if we are preserving the row-sum! + 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))) diff --git a/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 index 5b4f0a92..742a96fb 100644 --- a/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 @@ -110,6 +110,7 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& use amg_base_prec_type use amg_d_inner_mod, amg_protect_name => amg_daggrmat_smth_bld use amg_d_base_aggregator_mod +! use, intrinsic :: ieee_arithmetic implicit none @@ -192,7 +193,7 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (parms%aggr_filter == amg_filter_mat_) + filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_) ! ! naggr: number of local aggregates @@ -220,7 +221,7 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& ! \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)) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) end do !$OMP end parallel do end if @@ -257,8 +258,13 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& if (jd == -1) then ! if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i - else + else if (parms%aggr_filter == amg_filter_mat_) then + ! We perform filtering in the standard way assuming that A is an M-matrix acsrf%val(jd)=acsrf%val(jd)-tmp + else if (parms%aggr_filter == amg_filter_prow_mat_) then + ! We are probably doing l1-correction, hence we want to preserve the + ! row sum of the matrix: note the change in sign + acsrf%val(jd)=acsrf%val(jd)+tmp end if enddo !$OMP end parallel do @@ -266,7 +272,6 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& call acsrf%clean_zeros(info) end if - !$OMP parallel do private(i) schedule(static) do i=1,size(adiag) if (adiag(i) /= dzero) then @@ -278,11 +283,11 @@ subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - !if (do_l1correction) then - ! ! For l1-Jacobi this can be estimated with 1 - ! parms%aggr_omega_val = done - ! - if (parms%aggr_eig == amg_max_norm_) then + if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then + ! For l1-Jacobi this can be estimated with 1: + ! this makes sense only if we are preserving the row-sum! + 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))) diff --git a/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 index 1688b369..6caf61e5 100644 --- a/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 @@ -110,6 +110,7 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& use amg_base_prec_type use amg_s_inner_mod, amg_protect_name => amg_saggrmat_smth_bld use amg_s_base_aggregator_mod +! use, intrinsic :: ieee_arithmetic implicit none @@ -192,7 +193,7 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (parms%aggr_filter == amg_filter_mat_) + filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_) ! ! naggr: number of local aggregates @@ -220,7 +221,7 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& ! \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)) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) end do !$OMP end parallel do end if @@ -257,8 +258,13 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& if (jd == -1) then ! if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i - else + else if (parms%aggr_filter == amg_filter_mat_) then + ! We perform filtering in the standard way assuming that A is an M-matrix acsrf%val(jd)=acsrf%val(jd)-tmp + else if (parms%aggr_filter == amg_filter_prow_mat_) then + ! We are probably doing l1-correction, hence we want to preserve the + ! row sum of the matrix: note the change in sign + acsrf%val(jd)=acsrf%val(jd)+tmp end if enddo !$OMP end parallel do @@ -266,7 +272,6 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& call acsrf%clean_zeros(info) end if - !$OMP parallel do private(i) schedule(static) do i=1,size(adiag) if (adiag(i) /= szero) then @@ -278,11 +283,11 @@ subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - !if (do_l1correction) then - ! ! For l1-Jacobi this can be estimated with 1 - ! parms%aggr_omega_val = done - ! - if (parms%aggr_eig == amg_max_norm_) then + if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then + ! For l1-Jacobi this can be estimated with 1: + ! this makes sense only if we are preserving the row-sum! + 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))) diff --git a/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 index d35ac125..7d2cbe5f 100644 --- a/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 @@ -110,6 +110,7 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& use amg_base_prec_type use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_smth_bld use amg_z_base_aggregator_mod +! use, intrinsic :: ieee_arithmetic implicit none @@ -192,7 +193,7 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (parms%aggr_filter == amg_filter_mat_) + filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_) ! ! naggr: number of local aggregates @@ -220,7 +221,7 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& ! \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)) + adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i)) end do !$OMP end parallel do end if @@ -257,8 +258,13 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& if (jd == -1) then ! if (.not.do_l1correction) write(0,*) 'Wrong input: we need the diagonal!!!!', i - else + else if (parms%aggr_filter == amg_filter_mat_) then + ! We perform filtering in the standard way assuming that A is an M-matrix acsrf%val(jd)=acsrf%val(jd)-tmp + else if (parms%aggr_filter == amg_filter_prow_mat_) then + ! We are probably doing l1-correction, hence we want to preserve the + ! row sum of the matrix: note the change in sign + acsrf%val(jd)=acsrf%val(jd)+tmp end if enddo !$OMP end parallel do @@ -266,7 +272,6 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& call acsrf%clean_zeros(info) end if - !$OMP parallel do private(i) schedule(static) do i=1,size(adiag) if (adiag(i) /= zzero) then @@ -278,11 +283,11 @@ subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,& !$OMP end parallel do if (parms%aggr_omega_alg == amg_eig_est_) then - !if (do_l1correction) then - ! ! For l1-Jacobi this can be estimated with 1 - ! parms%aggr_omega_val = done - ! - if (parms%aggr_eig == amg_max_norm_) then + if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then + ! For l1-Jacobi this can be estimated with 1: + ! this makes sense only if we are preserving the row-sum! + 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)))