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

Merge branch 'l1-and-0-aggr' into repackage

Browse Source
This commit is contained in:
sfilippone
2024-11-12 10:07:44 +01:00
parent 053d5c4bc0 83ba79d7ae
commit d518f5eac7
37 changed files with 400 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files
amgprec/amg_base_prec_type.F90
+7 -4
View File
@@ -290,9 +290,10 @@ module amg_base_prec_type
!
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_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')
amgprec/amg_c_inner_mod.f90
+2 -1
View File
@@ -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(:)
amgprec/amg_d_inner_mod.f90
+2 -1
View File
@@ -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(:)
amgprec/amg_d_parmatch_aggregator_mod.F90
+4 -2
View File
@@ -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
amgprec/amg_s_inner_mod.f90
+2 -1
View File
@@ -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(:)
amgprec/amg_s_parmatch_aggregator_mod.F90
+4 -2
View File
@@ -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
amgprec/amg_z_inner_mod.f90
+2 -1
View File
@@ -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(:)
amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90
+8 -7
View File
@@ -177,13 +177,14 @@ 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_)
!!$
@@ -192,8 +193,8 @@ subroutine amg_c_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
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_
amgprec/impl/aggregator/amg_c_soc1_map_bld.F90
+2 -3
View File
@@ -250,7 +250,7 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
! we will not reset.
if (j>nr) cycle step1
if (ilaggr(j) > 0) cycle step1
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.czero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -357,7 +357,7 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.czero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -545,4 +545,3 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
return
end subroutine amg_c_soc1_map_bld
amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90
+11 -2
View File
@@ -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
amgprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90
+11 -3
View File
@@ -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()
amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90
+35 -5
View File
@@ -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,8 @@ 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,6 +278,10 @@ 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 (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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
@@ -259,7 +289,6 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
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 +351,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
amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90
+8 -7
View File
@@ -177,13 +177,14 @@ 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_)
!!$
@@ -192,8 +193,8 @@ subroutine amg_d_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
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_
amgprec/impl/aggregator/amg_d_parmatch_aggregator_mat_bld.F90
+10 -7
View File
@@ -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_
amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90
+29 -4
View File
@@ -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_) &
amgprec/impl/aggregator/amg_d_parmatch_unsmth_bld.F90
+10 -2
View File
@@ -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()
amgprec/impl/aggregator/amg_d_soc1_map_bld.F90
+2 -3
View File
@@ -250,7 +250,7 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
! we will not reset.
if (j>nr) cycle step1
if (ilaggr(j) > 0) cycle step1
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.dzero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -357,7 +357,7 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.dzero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -545,4 +545,3 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
return
end subroutine amg_d_soc1_map_bld
amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90
+11 -2
View File
@@ -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
amgprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90
+11 -3
View File
@@ -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()
amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90
+35 -5
View File
@@ -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,8 @@ 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,6 +278,10 @@ 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 (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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
@@ -259,7 +289,6 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
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 +351,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
amgprec/impl/aggregator/amg_s_dec_aggregator_mat_bld.f90
+8 -7
View File
@@ -177,13 +177,14 @@ 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_)
!!$
@@ -192,8 +193,8 @@ subroutine amg_s_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
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_
amgprec/impl/aggregator/amg_s_parmatch_aggregator_mat_bld.F90
+10 -7
View File
@@ -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_
amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90
+29 -4
View File
@@ -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_) &
amgprec/impl/aggregator/amg_s_parmatch_unsmth_bld.F90
+10 -2
View File
@@ -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()
amgprec/impl/aggregator/amg_s_soc1_map_bld.F90
+2 -3
View File
@@ -250,7 +250,7 @@ subroutine amg_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
! we will not reset.
if (j>nr) cycle step1
if (ilaggr(j) > 0) cycle step1
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.szero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -357,7 +357,7 @@ subroutine amg_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.szero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -545,4 +545,3 @@ subroutine amg_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
return
end subroutine amg_s_soc1_map_bld
amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90
+11 -2
View File
@@ -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
amgprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90
+11 -3
View File
@@ -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()
amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90
+35 -5
View File
@@ -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,8 @@ 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,6 +278,10 @@ 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 (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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
@@ -259,7 +289,6 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
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 +351,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
amgprec/impl/aggregator/amg_z_dec_aggregator_mat_bld.f90
+8 -7
View File
@@ -177,13 +177,14 @@ 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_)
!!$
@@ -192,8 +193,8 @@ subroutine amg_z_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
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_
amgprec/impl/aggregator/amg_z_soc1_map_bld.F90
+2 -3
View File
@@ -250,7 +250,7 @@ subroutine amg_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
! we will not reset.
if (j>nr) cycle step1
if (ilaggr(j) > 0) cycle step1
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.zzero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -357,7 +357,7 @@ subroutine amg_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
do k=1, nz
j = icol(k)
if ((1<=j).and.(j<=nr)) then
if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.zzero)) then
ip = ip + 1
icol(ip) = icol(k)
end if
@@ -545,4 +545,3 @@ subroutine amg_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
return
end subroutine amg_z_soc1_map_bld
amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90
+11 -2
View File
@@ -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
amgprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90
+11 -3
View File
@@ -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()
amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90
+35 -5
View File
@@ -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,8 @@ 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,6 +278,10 @@ 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 (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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
@@ -259,7 +289,6 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
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 +351,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
amgprec/impl/smoother/amg_d_poly_smoother_dmp.f90
+1 -3
View File
@@ -77,12 +77,10 @@ subroutine amg_d_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solve
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_poly',iam
write(fname(lname+1:lname+8),'(a,i3.3)') '_poly',iam
lname = lname + 8
! to be completed
! At base level do nothing for the smoother
if (allocated(sm%sv)) &
& call sm%sv%dump(desc,level,info,solver=solver,prefix=prefix,global_num=global_num)
amgprec/impl/smoother/amg_s_poly_smoother_dmp.f90
+1 -3
View File
@@ -77,12 +77,10 @@ subroutine amg_s_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solve
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_poly',iam
write(fname(lname+1:lname+8),'(a,i3.3)') '_poly',iam
lname = lname + 8
! to be completed
! At base level do nothing for the smoother
if (allocated(sm%sv)) &
& call sm%sv%dump(desc,level,info,solver=solver,prefix=prefix,global_num=global_num)
samples/advanced/pdegen/runs/amg_pde2d.inp
+1 -1
View File
@@ -20,7 +20,7 @@ FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r
1 ! degree for polynomial smoother
CHEB_4_OPT ! Polynomial variant
% Fields to be added for POLY
POLY_RHO_EST_POWER % POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER
POLY_RHO_EST_POWER ! POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER
% POLY_RHO_ESTIMATE_ITERATIONS default = 20
% POLY_RHO_BA set to value
1.0
samples/advanced/pdegen/runs/amg_pde3d.inp
+2 -2
View File
@@ -1,6 +1,6 @@
%%%%%%%%%%% General arguments % Lines starting with % are ignored.
CSR ! Storage format CSR COO JAD
0150 ! IDIM; domain size. Linear system size is IDIM**3
0050 ! IDIM; domain size. Linear system size is IDIM**3
POISSON ! PDECOEFF: POISSON, EXP, BOX, GAUSS
% ! Coefficients of the PDE
CG ! Iterative method:
@@ -20,7 +20,7 @@ POLY ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r
6 ! degree for polynomial smoother
CHEB_4_OPT ! Polynomial variant
% Fields to be added for POLY
POLY_RHO_EST_POWER % POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER
POLY_RHO_EST_POWER ! POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER
% POLY_RHO_ESTIMATE_ITERATIONS default = 20
% POLY_RHO_BA set to value
1.0