Added POLY smoothers, also in SAMPLES/ADVANCED

Poly-novrl
sfilippone 1 year ago
parent 737ebb9a96
commit 30a5c7be03

@ -21,7 +21,7 @@ DMODOBJS=amg_d_prec_type.o \
SMODOBJS=amg_s_prec_type.o amg_s_ilu_fact_mod.o \
amg_s_inner_mod.o amg_s_ilu_solver.o amg_s_diag_solver.o amg_s_jac_smoother.o amg_s_as_smoother.o \
amg_s_slu_solver.o amg_s_id_solver.o\
amg_s_poly_smoother.o amg_s_slu_solver.o amg_s_id_solver.o\
amg_s_base_solver_mod.o amg_s_base_smoother_mod.o amg_s_onelev_mod.o \
amg_s_gs_solver.o amg_s_mumps_solver.o amg_s_jac_solver.o \
amg_s_base_aggregator_mod.o \
@ -166,6 +166,7 @@ amg_dprecinit.o amg_dprecset.o: amg_d_diag_solver.o amg_d_ilu_solver.o \
amg_d_umf_solver.o amg_d_as_smoother.o amg_d_jac_smoother.o \
amg_d_id_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o
amg_d_poly_smoother.o: amg_d_base_smoother_mod.o amg_d_poly_coeff_mod.o
amg_s_poly_smoother.o: amg_s_base_smoother_mod.o amg_d_poly_coeff_mod.o
amg_s_mumps_solver.o amg_s_gs_solver.o amg_s_id_solver.o amg_s_slu_solver.o \
amg_s_diag_solver.o amg_s_ilu_solver.o amg_s_jac_solver.o: amg_s_base_solver_mod.o amg_s_prec_type.o

@ -687,10 +687,10 @@ contains
& ml_names(pm%ml_cycle)
select case (pm%ml_cycle)
case (amg_add_ml_)
write(iout,*) ' Number of smoother sweeps : ',&
write(iout,*) ' Number of smoother sweeps/degree : ',&
& pm%sweeps_pre
case (amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_, amg_kcycle_ml_, amg_kcyclesym_ml_)
write(iout,*) ' Number of smoother sweeps : pre: ',&
write(iout,*) ' Number of smoother sweeps/degree : pre: ',&
& pm%sweeps_pre ,' post: ', pm%sweeps_post
end select

@ -170,7 +170,7 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix)
call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix)
nswps = prec%precv(1)%parms%sweeps_pre
end if
write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps
write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps
write(iout_,*) trim(prefix_)
else if (nlev > 1) then

@ -98,6 +98,7 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
use amg_c_diag_solver
use amg_c_ilu_solver
use amg_c_gs_solver
#if defined(HAVE_SLU_)
use amg_c_slu_solver
#endif
@ -152,7 +153,6 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
if (info /= psb_success_) return
allocate(amg_c_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default()
case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
nlev_ = 1
ilev_ = 1

@ -170,7 +170,7 @@ subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity,prefix)
call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix)
nswps = prec%precv(1)%parms%sweeps_pre
end if
write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps
write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps
write(iout_,*) trim(prefix_)
else if (nlev > 1) then

@ -93,12 +93,13 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
use psb_base_mod
use amg_d_prec_mod, amg_protect_name => amg_dprecinit
use amg_d_jac_smoother
use amg_d_poly_smoother
use amg_d_as_smoother
use amg_d_id_solver
use amg_d_diag_solver
use amg_d_ilu_solver
use amg_d_gs_solver
use amg_d_poly_smoother
#if defined(HAVE_UMF_)
use amg_d_umf_solver
#endif
@ -156,7 +157,6 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
if (info /= psb_success_) return
allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default()
case ('POLY')
nlev_ = 1
ilev_ = 1
@ -165,7 +165,6 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
if (info /= psb_success_) return
allocate(amg_d_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default()
case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
nlev_ = 1
ilev_ = 1

@ -170,7 +170,7 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix)
call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix)
nswps = prec%precv(1)%parms%sweeps_pre
end if
write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps
write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps
write(iout_,*) trim(prefix_)
else if (nlev > 1) then

@ -98,6 +98,8 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
use amg_s_diag_solver
use amg_s_ilu_solver
use amg_s_gs_solver
use amg_s_poly_smoother
#if defined(HAVE_SLU_)
use amg_s_slu_solver
#endif
@ -152,7 +154,14 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
if (info /= psb_success_) return
allocate(amg_s_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default()
case ('POLY')
nlev_ = 1
ilev_ = 1
allocate(prec%precv(nlev_),stat=info)
allocate(amg_s_poly_smoother_type :: prec%precv(ilev_)%sm, stat=info)
if (info /= psb_success_) return
allocate(amg_s_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default()
case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
nlev_ = 1
ilev_ = 1

@ -170,7 +170,7 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix)
call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix)
nswps = prec%precv(1)%parms%sweeps_pre
end if
write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps
write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps
write(iout_,*) trim(prefix_)
else if (nlev > 1) then

@ -98,6 +98,7 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
use amg_z_diag_solver
use amg_z_ilu_solver
use amg_z_gs_solver
#if defined(HAVE_UMF_)
use amg_z_umf_solver
#endif
@ -155,7 +156,6 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
if (info /= psb_success_) return
allocate(amg_z_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default()
case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
nlev_ = 1
ilev_ = 1

@ -188,16 +188,11 @@ subroutine amg_c_base_onelev_csetc(lv,what,val,info,pos,idx)
case ('NONE','NOPREC','FACT_NONE')
call lv%set(amg_c_id_solver_mold,info,pos=pos)
case ('DIAG')
case ('DIAG','JACOBI')
call lv%set(amg_c_diag_solver_mold,info,pos=pos)
case ('JACOBI')
call lv%set(amg_c_jac_solver_mold,info,pos=pos)
case ('L1-DIAG')
case ('L1-DIAG','L1-JACOBI')
call lv%set(amg_c_l1_diag_solver_mold,info,pos=pos)
case ('L1-JACOBI')
call lv%set(amg_c_l1_jac_solver_mold,info,pos=pos)
case ('GS','FGS','FWGS')
call lv%set(amg_c_gs_solver_mold,info,pos=pos)

@ -43,9 +43,9 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
use amg_d_dec_aggregator_mod
use amg_d_symdec_aggregator_mod
use amg_d_parmatch_aggregator_mod
use amg_d_poly_smoother
use amg_d_jac_smoother
use amg_d_as_smoother
use amg_d_poly_smoother
use amg_d_diag_solver
use amg_d_l1_diag_solver
use amg_d_jac_solver
@ -85,7 +85,6 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
type(amg_d_jac_smoother_type) :: amg_d_jac_smoother_mold
type(amg_d_l1_jac_smoother_type) :: amg_d_l1_jac_smoother_mold
type(amg_d_as_smoother_type) :: amg_d_as_smoother_mold
type(amg_d_poly_smoother_type) :: amg_d_poly_smoother_mold
type(amg_d_diag_solver_type) :: amg_d_diag_solver_mold
type(amg_d_l1_diag_solver_type) :: amg_d_l1_diag_solver_mold
type(amg_d_jac_solver_type) :: amg_d_jac_solver_mold
@ -97,6 +96,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
type(amg_d_ainv_solver_type) :: amg_d_ainv_solver_mold
type(amg_d_invk_solver_type) :: amg_d_invk_solver_mold
type(amg_d_invt_solver_type) :: amg_d_invt_solver_mold
type(amg_d_poly_smoother_type) :: amg_d_poly_smoother_mold
#if defined(HAVE_UMF_)
type(amg_d_umf_solver_type) :: amg_d_umf_solver_mold
#endif
@ -161,7 +161,6 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
case ('POLY')
call lv%set(amg_d_poly_smoother_mold,info,pos=pos)
if (info == 0) call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos)
case ('GS','FWGS')
call lv%set(amg_d_jac_smoother_mold,info,pos='pre')
if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre')

@ -43,6 +43,7 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
use amg_s_dec_aggregator_mod
use amg_s_symdec_aggregator_mod
use amg_s_parmatch_aggregator_mod
use amg_s_poly_smoother
use amg_s_jac_smoother
use amg_s_as_smoother
use amg_s_diag_solver
@ -89,6 +90,7 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
type(amg_s_ainv_solver_type) :: amg_s_ainv_solver_mold
type(amg_s_invk_solver_type) :: amg_s_invk_solver_mold
type(amg_s_invt_solver_type) :: amg_s_invt_solver_mold
type(amg_s_poly_smoother_type) :: amg_s_poly_smoother_mold
#if defined(HAVE_SLU_)
type(amg_s_slu_solver_type) :: amg_s_slu_solver_mold
#endif
@ -144,6 +146,9 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
call lv%set(amg_s_as_smoother_mold,info,pos=pos)
if (info == 0) call lv%set(amg_s_ilu_solver_mold,info,pos=pos)
case ('POLY')
call lv%set(amg_s_poly_smoother_mold,info,pos=pos)
if (info == 0) call lv%set(amg_s_l1_diag_solver_mold,info,pos=pos)
case ('GS','FWGS')
call lv%set(amg_s_jac_smoother_mold,info,pos='pre')
if (info == 0) call lv%set(amg_s_gs_solver_mold,info,pos='pre')
@ -189,16 +194,11 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
case ('NONE','NOPREC','FACT_NONE')
call lv%set(amg_s_id_solver_mold,info,pos=pos)
case ('DIAG')
case ('DIAG','JACOBI')
call lv%set(amg_s_diag_solver_mold,info,pos=pos)
case ('JACOBI')
call lv%set(amg_s_jac_solver_mold,info,pos=pos)
case ('L1-DIAG')
case ('L1-DIAG','L1-JACOBI')
call lv%set(amg_s_l1_diag_solver_mold,info,pos=pos)
case ('L1-JACOBI')
call lv%set(amg_s_l1_jac_solver_mold,info,pos=pos)
case ('GS','FGS','FWGS')
call lv%set(amg_s_gs_solver_mold,info,pos=pos)

@ -200,16 +200,11 @@ subroutine amg_z_base_onelev_csetc(lv,what,val,info,pos,idx)
case ('NONE','NOPREC','FACT_NONE')
call lv%set(amg_z_id_solver_mold,info,pos=pos)
case ('DIAG')
case ('DIAG','JACOBI')
call lv%set(amg_z_diag_solver_mold,info,pos=pos)
case ('JACOBI')
call lv%set(amg_z_jac_solver_mold,info,pos=pos)
case ('L1-DIAG')
case ('L1-DIAG','L1-JACOBI')
call lv%set(amg_z_l1_diag_solver_mold,info,pos=pos)
case ('L1-JACOBI')
call lv%set(amg_z_l1_jac_solver_mold,info,pos=pos)
case ('GS','FGS','FWGS')
call lv%set(amg_z_gs_solver_mold,info,pos=pos)

@ -153,6 +153,17 @@ amg_s_jac_smoother_csetr.o \
amg_s_l1_jac_smoother_bld.o \
amg_s_l1_jac_smoother_descr.o \
amg_s_l1_jac_smoother_clone.o \
amg_s_poly_smoother_apply_vect.o \
amg_s_poly_smoother_bld.o \
amg_s_poly_smoother_cnv.o \
amg_s_poly_smoother_clone.o \
amg_s_poly_smoother_clone_settings.o \
amg_s_poly_smoother_clear_data.o \
amg_s_poly_smoother_descr.o \
amg_s_poly_smoother_dmp.o \
amg_s_poly_smoother_csetc.o \
amg_s_poly_smoother_cseti.o \
amg_s_poly_smoother_csetr.o \
amg_z_as_smoother_apply.o \
amg_z_as_smoother_apply_vect.o \
amg_z_as_smoother_bld.o \

@ -175,7 +175,7 @@ subroutine amg_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (info /= psb_success_) exit
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
call psb_geaxpby(cone,x,czero,r,r,desc_data,info)
call psb_geaxpby(cone,x,czero,r,desc_data,info)
call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info)
res = psb_genrm2(r,desc_data,info)
if( sm%printres ) then

@ -109,7 +109,7 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
if(.true..or.sm%checkres) then
if(sm%checkres) then
call psb_geall(r,desc_data,info)
call psb_geasb(r,desc_data,info)
resdenum = psb_genrm2(x,desc_data,info)
@ -159,10 +159,7 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
!!$ call psb_geaxpby(done,x,dzero,r,desc_data,info)
!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Jacobi smoother ',1,res
do i=1, sweeps-1
!
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
@ -176,10 +173,6 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
if (info /= psb_success_) exit
!!$ call psb_geaxpby(done,x,dzero,r,desc_data,info)
!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Jacobi smoother ',i+1,res
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
call psb_geaxpby(done,x,dzero,r,desc_data,info)

@ -120,7 +120,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& a_err='invalid wv size in smoother_apply')
goto 9999
end if
sm%pdegree = sweeps
associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4))
call psb_geaxpby(done,x,dzero,r,desc_data,info)
@ -135,7 +135,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! b == x
! x == tx
!
do i=1, sm%pdegree
do i=1, sweeps
! B r_{k-1}
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
cz = (2*i*done-3)/(2*i*done+done)
@ -163,7 +163,15 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! b == x
! x == tx
!
do i=1, sm%pdegree
if (allocated(sm%poly_beta)) then
if (size(sm%poly_beta) /= sweeps) deallocate(sm%poly_beta)
end if
if (.not.allocated(sm%poly_beta)) then
call psb_realloc(sweeps,sm%poly_beta,info)
sm%poly_beta(1:sweeps) = amg_d_poly_beta_mat(1:sweeps,sweeps)
end if
do i=1, sweeps
! B r_{k-1}
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
cz = (2*i*done-3)/(2*i*done+done)
@ -190,6 +198,8 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! b == x
! x == tx
!
sm%cf_a = amg_d_poly_a_vect(sweeps)
theta = (done+sm%cf_a)/2
delta = (done-sm%cf_a)/2
sigma = theta/delta
@ -198,7 +208,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
call psb_geaxpby((done/sm%rho_ba),ty,dzero,r,desc_data,info)
call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info)
! tz == d
do i=1, sm%pdegree
do i=1, sweeps
! x_{k+1} = x_k + d_k
call psb_geaxpby(done,tz,done,tx,desc_data,info)
!

@ -74,7 +74,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
n_col = desc_a%get_local_cols()
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if (.false.) then
select case(sm%variant)
case(amg_poly_lottes_)
! do nothing
@ -106,7 +106,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='invalid sm%variant')
goto 9999
end select
end if
sm%pa => a
if (.not.allocated(sm%sv)) then
info = psb_err_internal_error_

@ -44,7 +44,7 @@ subroutine amg_d_poly_smoother_clear_data(sm,info)
class(amg_d_poly_smoother_type), intent(inout) :: sm
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act
character(len=20) :: name='d_poly_smoother_clear_data'
character(len=20) :: name='amg_d_poly_smoother_clear_data'
call psb_erractionsave(err_act)

@ -56,6 +56,12 @@ subroutine amg_d_poly_smoother_clone_settings(sm,smout,info)
smout%pa => null()
smout%pdegree = sm%pdegree
smout%variant = sm%variant
smout%cf_a = sm%cf_a
smout%rho_ba = sm%rho_ba
smout%rho_estimate = sm%rho_estimate
smout%rho_estimate_iterations = sm%rho_estimate_iterations
smout%poly_beta = sm%poly_beta
if (allocated(smout%sv)) then
if (.not.same_type_as(sm%sv,smout%sv)) then

@ -38,7 +38,7 @@
subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx)
use psb_base_mod
use amg_d_poly_smoother, amg_protect_nam => amg_d_poly_smoother_csetc
use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_csetc
Implicit None
! Arguments
class(amg_d_poly_smoother_type), intent(inout) :: sm

@ -40,8 +40,6 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix)
use psb_base_mod
use amg_d_diag_solver
use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_descr
use amg_d_diag_solver
use amg_d_gs_solver
Implicit None
@ -81,18 +79,18 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix)
select case(sm%variant)
case(amg_poly_lottes_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES'
write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
!write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
case(amg_poly_lottes_beta_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES_BETA'
write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
!write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), ' Coefficients: ',sm%poly_beta(1:sm%pdegree)
!if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), ' Coefficients: ',sm%poly_beta(1:sm%pdegree)
case(amg_poly_new_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_NEW'
write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
!write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a
!write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a
case default
write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???'
end select

@ -175,7 +175,7 @@ subroutine amg_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (info /= psb_success_) exit
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
call psb_geaxpby(sone,x,szero,r,desc_data,info)
call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
res = psb_genrm2(r,desc_data,info)
if( sm%printres ) then

@ -0,0 +1,453 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,wv,info,init,initu)
use psb_base_mod
use amg_s_diag_solver
use psb_base_krylov_conv_mod, only : log_conv
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_apply_vect
implicit none
type(psb_desc_type), intent(in) :: desc_data
class(amg_s_poly_smoother_type), intent(inout) :: sm
type(psb_s_vect_type),intent(inout) :: x
type(psb_s_vect_type),intent(inout) :: y
real(psb_spk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps
real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu
!
integer(psb_ipk_) :: n_row,n_col
type(psb_s_vect_type) :: tx, ty, tz, r
real(psb_spk_), pointer :: aux(:)
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me, i, err_act
character :: trans_, init_
real(psb_spk_) :: res, resdenum
character(len=20) :: name='d_poly_smoother_apply_v'
call psb_erractionsave(err_act)
info = psb_success_
ctxt = desc_data%get_context()
call psb_info(ctxt,me,np)
if (present(init)) then
init_ = psb_toupper(init)
else
init_='Z'
end if
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols()
if (4*n_col <= size(work)) then
aux => work(:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,&
& i_err=(/4*n_col,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
endif
!!$ if (me == 0) write(0,*) name,' Unimplemented apply_vect '
!!$ info =psb_err_internal_error_
!!$ call psb_errpush(info,&
!!$ & name,a_err='Error in sub_aply Polynomial not implemented')
!!$ goto 9999
if (size(wv) < 4) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size in smoother_apply')
goto 9999
end if
sm%pdegree = sweeps
associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4))
call psb_geaxpby(sone,x,szero,r,desc_data,info)
call tx%zero()
call ty%zero()
call tz%zero()
select case(sm%variant)
case(amg_poly_lottes_)
block
real(psb_spk_) :: cz, cr
! b == x
! x == tx
!
do i=1, sweeps
! B r_{k-1}
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
cz = (2*i*sone-3)/(2*i*sone+sone)
cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
! r_k = b-Ax_k = x -A tx
call psb_geaxpby(sone,tz,sone,tx,desc_data,info)
if (.false.) then
call psb_geaxpby(sone,x,szero,r,desc_data,info)
call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_)
else
call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
end if
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Polynomial smoother ',i,res
! x_k = x_{k-1} + z_k
end do
end block
case(amg_poly_lottes_beta_)
block
real(psb_spk_) :: cz, cr
! b == x
! x == tx
!
if (allocated(sm%poly_beta)) then
if (size(sm%poly_beta) /= sweeps) deallocate(sm%poly_beta)
end if
if (.not.allocated(sm%poly_beta)) then
call psb_realloc(sweeps,sm%poly_beta,info)
sm%poly_beta(1:sweeps) = amg_d_poly_beta_mat(1:sweeps,sweeps)
end if
do i=1, sweeps
! B r_{k-1}
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
cz = (2*i*sone-3)/(2*i*sone+sone)
cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
! r_k = b-Ax_k = x -A tx
call psb_geaxpby(sm%poly_beta(i),tz,sone,tx,desc_data,info)
if (.false.) then
call psb_geaxpby(sone,x,szero,r,desc_data,info)
call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_)
else
call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
end if
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Polynomial smoother ',i,res
! x_k = x_{k-1} + z_k
end do
end block
case(amg_poly_new_)
block
real(psb_spk_) :: sigma, theta, delta, rho_old, rho
! b == x
! x == tx
!
sm%cf_a = amg_d_poly_a_vect(sweeps)
theta = (sone+sm%cf_a)/2
delta = (sone-sm%cf_a)/2
sigma = theta/delta
rho_old = sone/sigma
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info)
call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info)
! tz == d
do i=1, sweeps
! x_{k+1} = x_k + d_k
call psb_geaxpby(sone,tz,sone,tx,desc_data,info)
!
! r_{k-1} = r_k - (1/rho(BA)) B A d_k
call psb_spmm(sone,sm%pa,tz,szero,ty,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(-sone,ty,sone,r,desc_data,trans_,aux,wv(5:),info,init='Z')
!
! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
rho = sone/(2*sigma - rho_old)
call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info)
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Polynomial smoother ',i,res
! x_k = x_{k-1} + z_k
end do
end block
case default
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='wrong polynomial variant')
goto 9999
end select
if (info == psb_success_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='polynomial smoother')
goto 9999
end if
end associate
!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
!!$ ! if .not.sv%is_iterative, there's no need to pass init
!!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)
!!$
!!$ if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,&
!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999
!!$ endif
!!$
!!$ else if (sweeps >= 0) then
!!$ select type (smsv => sm%sv)
!!$ class is (amg_s_diag_solver_type)
!!$ !
!!$ ! This means we are dealing with a pure Jacobi smoother/solver.
!!$ !
!!$ associate(tx => wv(1), ty => wv(2))
!!$ select case (init_)
!!$ case('Z')
!!$
!!$ call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
!!$
!!$ case('Y')
!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info)
!!$ call psb_geaxpby(sone,y,szero,ty,desc_data,info)
!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
!!$
!!$ case('U')
!!$ if (.not.present(initu)) then
!!$ call psb_errpush(psb_err_internal_error_,name,&
!!$ & a_err='missing initu to smoother_apply')
!!$ goto 9999
!!$ end if
!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info)
!!$ call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
!!$
!!$ case default
!!$ call psb_errpush(psb_err_internal_error_,name,&
!!$ & a_err='wrong init to smoother_apply')
!!$ goto 9999
!!$ end select
!!$
!!$ do i=1, sweeps-1
!!$ !
!!$ ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
!!$ ! where is the diagonal and A the matrix.
!!$ !
!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info)
!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
!!$
!!$ if (info /= psb_success_) exit
!!$
!!$ call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
!!$
!!$ if (info /= psb_success_) exit
!!$
!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
!!$ call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
!!$ call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
!!$ res = psb_genrm2(r,desc_data,info)
!!$ if( sm%printres ) then
!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
!!$ end if
!!$ if ( res < sm%tol*resdenum ) then
!!$ if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) &
!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
!!$ exit
!!$ end if
!!$ end if
!!$
!!$ end do
!!$
!!$
!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
!!$
!!$ if (info /= psb_success_) then
!!$ info=psb_err_internal_error_
!!$ call psb_errpush(info,name,&
!!$ & a_err='subsolve with Jacobi sweeps > 1')
!!$ goto 9999
!!$ end if
!!$
!!$
!!$ end associate
!!$
!!$ class default
!!$ !
!!$ !
!!$ ! Apply multiple sweeps of a block-Jacobi solver
!!$ ! to compute an approximate solution of a linear system.
!!$ !
!!$ !
!!$ if (size(wv) < 2) then
!!$ info = psb_err_internal_error_
!!$ call psb_errpush(info,name,&
!!$ & a_err='invalid wv size in smoother_apply')
!!$ goto 9999
!!$ end if
!!$ associate(tx => wv(1), ty => wv(2))
!!$
!!$ !
!!$ ! Unroll the first iteration and fold it inside SELECT CASE
!!$ ! this will save one AXPBY and one SPMM when INIT=Z, and will be
!!$ ! significant when sweeps=1 (a common case)
!!$ !
!!$ select case (init_)
!!$ case('Z')
!!$
!!$ call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z')
!!$
!!$ case('Y')
!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info)
!!$ call psb_geaxpby(sone,y,szero,ty,desc_data,info)
!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
!!$
!!$ case('U')
!!$ if (.not.present(initu)) then
!!$ call psb_errpush(psb_err_internal_error_,name,&
!!$ & a_err='missing initu to smoother_apply')
!!$ goto 9999
!!$ end if
!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info)
!!$ call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
!!$
!!$ case default
!!$ call psb_errpush(psb_err_internal_error_,name,&
!!$ & a_err='wrong init to smoother_apply')
!!$ goto 9999
!!$ end select
!!$
!!$ do i=1, sweeps-1
!!$ !
!!$ ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
!!$ ! block diagonal part and the remaining part of the local matrix
!!$ ! and Y(j) is the approximate solution at sweep j.
!!$ !
!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info)
!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
!!$
!!$ if (info /= psb_success_) exit
!!$
!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
!!$
!!$ if (info /= psb_success_) exit
!!$
!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
!!$ call psb_geaxpby(sone,x,szero,r,r,desc_data,info)
!!$ call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info)
!!$ res = psb_genrm2(r,desc_data,info)
!!$ if( sm%printres ) then
!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol)
!!$ end if
!!$ if (res < sm%tol*resdenum ) then
!!$ if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) &
!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol)
!!$ exit
!!$ end if
!!$ end if
!!$
!!$ end do
!!$
!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
!!$
!!$ if (info /= psb_success_) then
!!$ info=psb_err_internal_error_
!!$ call psb_errpush(info,name,&
!!$ & a_err='subsolve with Jacobi sweeps > 1')
!!$ goto 9999
!!$ end if
!!$
!!$ end associate
!!$ end select
!!$
!!$ else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
!!$ goto 9999
!!$
!!$ endif
!!$
if (.not.(4*n_col <= size(work))) then
deallocate(aux)
endif
!!$ if(sm%checkres) then
!!$ call psb_gefree(r,desc_data,info)
!!$ end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_apply_vect

@ -0,0 +1,180 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
use psb_base_mod
use amg_s_diag_solver
use amg_s_l1_diag_solver
use amg_d_poly_coeff_mod
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_bld
Implicit None
! Arguments
type(psb_sspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(inout) :: desc_a
class(amg_s_poly_smoother_type), intent(inout) :: sm
integer(psb_ipk_), intent(out) :: info
class(psb_s_base_sparse_mat), intent(in), optional :: amold
class(psb_s_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
type(psb_sspmat_type) :: tmpa
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros
type(psb_ctxt_type) :: ctxt
real(psb_spk_), allocatable :: da(:), dsv(:)
integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level
character(len=20) :: name='d_poly_smoother_bld', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
n_row = desc_a%get_local_rows()
n_col = desc_a%get_local_cols()
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if (.false.) then
select case(sm%variant)
case(amg_poly_lottes_)
! do nothing
case(amg_poly_lottes_beta_)
if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
call psb_realloc(sm%pdegree,sm%poly_beta,info)
sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
else
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid sm%degree for poly_beta')
goto 9999
end if
case(amg_poly_new_)
if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
!Ok
sm%cf_a = amg_d_poly_a_vect(sm%pdegree)
else
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid sm%degree for poly_a')
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid sm%variant')
goto 9999
end select
end if
sm%pa => a
if (.not.allocated(sm%sv)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='unallocated sm%sv')
goto 9999
end if
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='sv%build')
goto 9999
end if
!!$ if (.false.) then
!!$ select type(ssv => sm%sv)
!!$ class is(amg_s_l1_diag_solver_type)
!!$ da = a%arwsum(info)
!!$ dsv = ssv%dv%get_vect()
!!$ sm%rho_ba = maxval(da(1:n_row)*dsv(1:n_row))
!!$ class default
!!$ write(0,*) 'PolySmoother BUILD: only L1-Jacobi/L1-DIAG for now ',ssv%get_fmt()
!!$ sm%rho_ba = sone
!!$ end select
!!$ else
if (sm%rho_ba <= szero) then
select case(sm%rho_estimate)
case(amg_poly_rho_est_power_)
block
type(psb_s_vect_type) :: tq, tt, tz,wv(2)
real(psb_spk_) :: znrm, lambda
real(psb_spk_),allocatable :: work(:)
integer(psb_ipk_) :: i, n_cols
n_cols = desc_a%get_local_cols()
allocate(work(4*n_cols))
call psb_geasb(tz,desc_a,info,mold=vmold,scratch=.true.)
call psb_geasb(tt,desc_a,info,mold=vmold,scratch=.true.)
call psb_geasb(wv(1),desc_a,info,mold=vmold,scratch=.true.)
call psb_geasb(wv(2),desc_a,info,mold=vmold,scratch=.true.)
call psb_geall(tq,desc_a,info)
call tq%set(sone)
call psb_geasb(tq,desc_a,info,mold=vmold)
call psb_spmm(sone,a,tq,szero,tt,desc_a,info) !
call sm%sv%apply_v(sone,tt,szero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = BA q_k
do i=1,sm%rho_estimate_iterations
znrm = psb_genrm2(tz,desc_a,info) ! znrm = |z_k|_2
call psb_geaxpby((sone/znrm),tz,szero,tq,desc_a,info) ! q_k = z_k/znrm
call psb_spmm(sone,a,tq,szero,tt,desc_a,info) ! t_{k+1} = BA q_k
call sm%sv%apply_v(sone,tt,szero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = B t_{k+1}
lambda = psb_gedot(tq,tz,desc_a,info) ! lambda = q_k^T z_{k+1} = q_k^T BA q_k
!write(0,*) 'BLD: lambda estimate ',i,lambda
end do
sm%rho_ba = lambda
end block
case default
write(0,*) ' Unknown algorithm for RHO(BA) estimate, defaulting to a value of 1.0 '
sm%rho_ba = sone
end select
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_bld

@ -0,0 +1,70 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_clear_data(sm,info)
use psb_base_mod
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_clear_data
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(inout) :: sm
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act
character(len=20) :: name='amg_s_poly_smoother_clear_data'
call psb_erractionsave(err_act)
info = 0
sm%pdegree = 0
if (allocated(sm%poly_beta)) deallocate(sm%poly_beta)
sm%pa => null()
if ((info==0).and.allocated(sm%sv)) then
call sm%sv%clear_data(info)
end if
if (info /= 0) then
info = psb_err_internal_error_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_clear_data

@ -0,0 +1,90 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_clone(sm,smout,info)
use psb_base_mod
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_clone
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(inout) :: sm
class(amg_s_base_smoother_type), allocatable, intent(inout) :: smout
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
info=psb_success_
call psb_erractionsave(err_act)
if (allocated(smout)) then
call smout%free(info)
if (info == psb_success_) deallocate(smout, stat=info)
end if
if (info == psb_success_) &
& allocate(amg_s_poly_smoother_type :: smout, stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
goto 9999
end if
select type(smo => smout)
type is (amg_s_poly_smoother_type)
smo%pdegree = sm%pdegree
smo%rho_ba = sm%rho_ba
smo%poly_beta = sm%poly_beta
smo%pa => sm%pa
if ((info==psb_success_).and.(allocated(sm%sv))) then
allocate(smout%sv,mold=sm%sv,stat=info)
if (info == psb_success_) call sm%sv%clone(smo%sv,info)
end if
class default
info = psb_err_internal_error_
end select
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_clone

@ -0,0 +1,102 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! asd on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_clone_settings(sm,smout,info)
use psb_base_mod
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_clone_settings
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(inout) :: sm
class(amg_s_base_smoother_type), intent(inout) :: smout
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: err_act
character(len=20) :: name='d_poly_smoother_clone_settings'
call psb_erractionsave(err_act)
info = psb_success_
select type(smout)
class is(amg_s_poly_smoother_type)
smout%pa => null()
smout%pdegree = sm%pdegree
smout%variant = sm%variant
smout%cf_a = sm%cf_a
smout%rho_ba = sm%rho_ba
smout%rho_estimate = sm%rho_estimate
smout%rho_estimate_iterations = sm%rho_estimate_iterations
smout%poly_beta = sm%poly_beta
if (allocated(smout%sv)) then
if (.not.same_type_as(sm%sv,smout%sv)) then
call smout%sv%free(info)
if (info == 0) deallocate(smout%sv,stat=info)
end if
end if
if (info /= 0) then
info = psb_err_internal_error_
else
if (allocated(smout%sv)) then
if (same_type_as(sm%sv,smout%sv)) then
call sm%sv%clone_settings(smout%sv,info)
else
info = psb_err_internal_error_
end if
else
allocate(smout%sv,mold=sm%sv,stat=info)
if (info == 0) call sm%sv%clone_settings(smout%sv,info)
if (info /= 0) info = psb_err_internal_error_
end if
end if
class default
info = psb_err_internal_error_
end select
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_clone_settings

@ -0,0 +1,77 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_cnv(sm,info,amold,vmold,imold)
use psb_base_mod
use amg_s_diag_solver
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_cnv
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(inout) :: sm
integer(psb_ipk_), intent(out) :: info
class(psb_s_base_sparse_mat), intent(in), optional :: amold
class(psb_s_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
integer(psb_ipk_) :: i, err_act, debug_unit, debug_level
character(len=20) :: name='d_poly_smoother_cnv', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
if (allocated(sm%sv)) &
& call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='solver cnv')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_cnv

@ -0,0 +1,76 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_csetc(sm,what,val,info,idx)
use psb_base_mod
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_csetc
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(inout) :: sm
character(len=*), intent(in) :: what
character(len=*), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: idx
integer(psb_ipk_) :: err_act, ival
character(len=20) :: name='d_poly_smoother_csetc'
info = psb_success_
call psb_erractionsave(err_act)
select case(psb_toupper(trim(what)))
case('POLY_VARIANT')
call sm%set(what,amg_stringval(val),info,idx=idx)
case('POLY_RHO_ESTIMATE')
call sm%set(what,amg_stringval(val),info,idx=idx)
case default
call sm%amg_s_base_smoother_type%set(what,val,info,idx=idx)
end select
if (info /= psb_success_) then
info = psb_err_from_subroutine_
call psb_errpush(info, name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_csetc

@ -0,0 +1,92 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_cseti(sm,what,val,info,idx)
use psb_base_mod
use amg_s_poly_smoother, amg_protect_nam => amg_s_poly_smoother_cseti
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(inout) :: sm
character(len=*), intent(in) :: what
integer(psb_ipk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: idx
integer(psb_ipk_) :: err_act
character(len=20) :: name='d_poly_smoother_cseti'
info = psb_success_
call psb_erractionsave(err_act)
select case(psb_toupper(what))
case('POLY_DEGREE')
sm%pdegree = val
case('POLY_VARIANT')
select case(val)
case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_)
sm%variant = val
case default
write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val
sm%variant = amg_poly_lottes_
end select
case('POLY_RHO_ESTIMATE')
select case(val)
case (amg_poly_rho_est_power_)
sm%rho_estimate = val
case default
write(0,*) 'Invalid choice for POLY_RHO_ESTIMATE, defaulting to amg_poly_rho_power'
sm%variant = amg_poly_rho_est_power_
end select
case('POLY_RHO_ESTIMATE_ITERATIONS')
if (val>0) then
sm%rho_estimate_iterations = val
else
write(0,*) 'Invalid choice for POLY_RHO_ESTIMATE_ITERATIONS, defaulting to 20'
sm%variant = 20
end if
case default
call sm%amg_s_base_smoother_type%set(what,val,info,idx=idx)
end select
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_cseti

@ -0,0 +1,74 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Daniela di Serafino
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_csetr(sm,what,val,info,idx)
use psb_base_mod
use amg_s_poly_smoother, amg_protect_nam => amg_s_poly_smoother_csetr
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(inout) :: sm
character(len=*), intent(in) :: what
real(psb_spk_), intent(in) :: val
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: idx
integer(psb_ipk_) :: err_act
character(len=20) :: name='d_poly_smoother_csetr'
info = psb_success_
call psb_erractionsave(err_act)
select case(psb_toupper(what))
case('POLY_RHO_BA')
if ((szero<val).and.(val<=sone)) then
sm%rho_ba = val
else
write(0,*) 'Invalid choice for POLY_RHO_BA, defaulting to compute estimate'
sm%rho_ba = -sone
end if
case default
call sm%amg_s_base_smoother_type%set(what,val,info,idx=idx)
end select
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_csetr

@ -0,0 +1,108 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_descr(sm,info,iout,coarse,prefix)
use psb_base_mod
use amg_s_diag_solver
use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_descr
Implicit None
! Arguments
class(amg_s_poly_smoother_type), intent(in) :: sm
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: iout
logical, intent(in), optional :: coarse
character(len=*), intent(in), optional :: prefix
! Local variables
integer(psb_ipk_) :: err_act
character(len=20), parameter :: name='amg_s_poly_smoother_descr'
integer(psb_ipk_) :: iout_
logical :: coarse_
character(1024) :: prefix_
call psb_erractionsave(err_act)
info = psb_success_
if (present(coarse)) then
coarse_ = coarse
else
coarse_ = .false.
end if
if (present(iout)) then
iout_ = iout
else
iout_ = psb_out_unit
endif
if (present(prefix)) then
prefix_ = prefix
else
prefix_ = ""
end if
write(iout_,*) trim(prefix_), ' Polynomial smoother '
select case(sm%variant)
case(amg_poly_lottes_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES'
!write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
case(amg_poly_lottes_beta_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES_BETA'
!write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
!if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), ' Coefficients: ',sm%poly_beta(1:sm%pdegree)
case(amg_poly_new_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_NEW'
!write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree
write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
!write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a
case default
write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???'
end select
if (allocated(sm%sv)) then
write(iout_,*) trim(prefix_), ' Local solver details:'
call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix)
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine amg_s_poly_smoother_descr

@ -0,0 +1,90 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
subroutine amg_s_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solver,global_num)
use psb_base_mod
use amg_s_poly_smoother, amg_protect_nam => amg_s_poly_smoother_dmp
implicit none
class(amg_s_poly_smoother_type), intent(in) :: sm
type(psb_desc_type), intent(in) :: desc
integer(psb_ipk_), intent(in) :: level
integer(psb_ipk_), intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver, global_num
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
integer(psb_lpk_), allocatable :: iv(:)
logical :: smoother_, global_num_
! len of prefix_
info = 0
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_smth_d"
end if
ctxt = desc%get_context()
call psb_info(ctxt,iam,np)
if (present(smoother)) then
smoother_ = smoother
else
smoother_ = .false.
end if
if (present(global_num)) then
global_num_ = global_num
else
global_num_ = .false.
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(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)
end subroutine amg_s_poly_smoother_dmp

@ -175,7 +175,7 @@ subroutine amg_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (info /= psb_success_) exit
if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
call psb_geaxpby(zone,x,zzero,r,r,desc_data,info)
call psb_geaxpby(zone,x,zzero,r,desc_data,info)
call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info)
res = psb_genrm2(r,desc_data,info)
if( sm%printres ) then

@ -8,23 +8,27 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUD
LINKOPT=
EXEDIR=./runs
DGEN2D=amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o
DGEN3D=amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o
SGEN2D=amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o
SGEN3D=amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o amg_s_pde3d_box_mod.o
all: amg_s_pde3d amg_d_pde3d amg_s_pde2d amg_d_pde2d
amg_d_pde3d: amg_d_pde3d.o amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o data_input.o
$(FLINK) $(LINKOPT) amg_d_pde3d.o amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
amg_d_pde3d: amg_d_pde3d.o amg_d_genpde_mod.o $(DGEN3D) data_input.o
$(FLINK) $(LINKOPT) amg_d_pde3d.o amg_d_genpde_mod.o $(DGEN3D) data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
/bin/mv amg_d_pde3d $(EXEDIR)
amg_s_pde3d: amg_s_pde3d.o amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o data_input.o
$(FLINK) $(LINKOPT) amg_s_pde3d.o amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
amg_s_pde3d: amg_s_pde3d.o amg_s_genpde_mod.o $(SGEN3D) data_input.o
$(FLINK) $(LINKOPT) amg_s_pde3d.o amg_s_genpde_mod.o $(SGEN3D) data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
/bin/mv amg_s_pde3d $(EXEDIR)
amg_d_pde2d: amg_d_pde2d.o amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o data_input.o
$(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
amg_d_pde2d: amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o
$(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
/bin/mv amg_d_pde2d $(EXEDIR)
amg_s_pde2d: amg_s_pde2d.o amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o data_input.o
$(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
amg_s_pde2d: amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o
$(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
/bin/mv amg_s_pde2d $(EXEDIR)
amg_d_pde3d_rebld: amg_d_pde3d_rebld.o data_input.o
@ -33,10 +37,10 @@ amg_d_pde3d_rebld: amg_d_pde3d_rebld.o data_input.o
amg_d_pde3d.o amg_s_pde3d.o amg_d_pde2d.o amg_s_pde2d.o: data_input.o
amg_d_pde3d.o: amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o
amg_s_pde3d.o: amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o
amg_d_pde2d.o: amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o
amg_s_pde2d.o: amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o
amg_d_pde3d.o: amg_d_genpde_mod.o $(DGEN3D)
amg_s_pde3d.o: amg_s_genpde_mod.o $(SGEN3D)
amg_d_pde2d.o: amg_d_genpde_mod.o $(DGEN2D)
amg_s_pde2d.o: amg_s_genpde_mod.o $(SGEN2D)
check: all
cd runs && ./amg_d_pde2d <amg_pde2d.inp && ./amg_s_pde2d<amg_pde2d.inp

@ -72,6 +72,7 @@ program amg_d_pde2d
use amg_d_pde2d_base_mod
use amg_d_pde2d_exp_mod
use amg_d_pde2d_box_mod
use amg_d_pde2d_gauss_mod
use amg_d_genpde_mod
#if defined(OPENMP)
use omp_lib
@ -124,16 +125,16 @@ program amg_d_pde2d
integer(psb_ipk_) :: outer_sweeps ! number of outer sweeps: sweeps for 1-level,
! AMG cycles for ML
! general AMG data
character(len=16) :: mlcycle ! AMG cycle type
character(len=32) :: mlcycle ! AMG cycle type
integer(psb_ipk_) :: maxlevs ! maximum number of levels in AMG preconditioner
! AMG aggregation
character(len=16) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED
character(len=16) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC
character(len=16) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP
character(len=32) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED
character(len=32) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC
character(len=32) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP
integer(psb_ipk_) :: aggr_size ! Requested size of the aggregates for MATCHBOXP
character(len=16) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE
character(len=16) :: aggr_filter ! filtering: FILTER, NO_FILTER
character(len=32) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE
character(len=32) :: aggr_filter ! filtering: FILTER, NO_FILTER
real(psb_dpk_) :: mncrratio ! minimum aggregation ratio
real(psb_dpk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector
integer(psb_ipk_) :: thrvsz ! size of threshold vector
@ -141,39 +142,43 @@ program amg_d_pde2d
integer(psb_ipk_) :: csizepp ! minimum size of coarsest matrix per process
! AMG smoother or pre-smoother; also 1-lev preconditioner
character(len=16) :: smther ! (pre-)smoother type: BJAC, AS
character(len=32) :: smther ! (pre-)smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps
integer(psb_ipk_) :: degree ! degree for polynomial smoother
character(len=32) :: pvariant ! polynomial variant
integer(psb_ipk_) :: novr ! number of overlap layers
character(len=16) :: restr ! restriction over application of AS
character(len=16) :: prol ! prolongation over application of AS
character(len=16) :: solve ! local subsolver type: ILU, MILU, ILUT,
character(len=32) :: restr ! restriction over application of AS
character(len=32) :: prol ! prolongation over application of AS
character(len=32) :: solve ! local subsolver type: ILU, MILU, ILUT,
! UMF, MUMPS, SLU, FWGS, BWGS, JAC
integer(psb_ipk_) :: ssweeps ! inner solver sweeps
character(len=16) :: variant ! AINV variant: LLK, etc
character(len=32) :: variant ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill ! Inverse fill-in for INVK
real(psb_dpk_) :: thr ! threshold for ILUT factorization
! AMG post-smoother; ignored by 1-lev preconditioner
character(len=16) :: smther2 ! post-smoother type: BJAC, AS
character(len=32) :: smther2 ! post-smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps
integer(psb_ipk_) :: degree2 ! degree for polynomial smoother
character(len=32) :: pvariant2 ! polynomial variant
integer(psb_ipk_) :: novr2 ! number of overlap layers
character(len=16) :: restr2 ! restriction over application of AS
character(len=16) :: prol2 ! prolongation over application of AS
character(len=16) :: solve2 ! local subsolver type: ILU, MILU, ILUT,
character(len=32) :: restr2 ! restriction over application of AS
character(len=32) :: prol2 ! prolongation over application of AS
character(len=32) :: solve2 ! local subsolver type: ILU, MILU, ILUT,
! UMF, MUMPS, SLU, FWGS, BWGS, JAC
integer(psb_ipk_) :: ssweeps2 ! inner solver sweeps
character(len=16) :: variant2 ! AINV variant: LLK, etc
character(len=32) :: variant2 ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill2 ! Inverse fill-in for INVK
real(psb_dpk_) :: thr2 ! threshold for ILUT factorization
! coarsest-level solver
character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST
character(len=16) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr.
character(len=32) :: cmat ! coarsest matrix layout: REPL, DIST
character(len=32) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr.
! mat.); UMF, MUMPS, SLU, ILU, ILUT, MILU
! (repl. mat.)
character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT,
character(len=32) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT,
! MILU, UMF, MUMPS, SLU
integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization
real(psb_dpk_) :: cthres ! threshold for ILUT factorization
@ -241,13 +246,16 @@ program amg_d_pde2d
select case(psb_toupper(trim(pdecoeff)))
case("CONST")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info)
& a1_base,a2_base,b1_base,b2_base,c_base,g_base,info)
case("EXP")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_exp,a2_exp,b1_exp,b2_exp,c_exp,g_exp,info)
case("BOX")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_box,a2_box,b1_box,b2_box,c_box,g_box,info)
case("GAUSS")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_gauss,a2_gauss,b1_gauss,b2_gauss,c_gauss,g_gauss,info)
case default
info=psb_err_from_subroutine_
ch_err='amg_gen_pdecoeff'
@ -281,10 +289,12 @@ program amg_d_pde2d
! 1-level sweeps from "outer_sweeps"
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
case ('BJAC')
case ('BJAC','POLY')
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('sub_solve', p_choice%solve, info)
call prec%set('solver_sweeps', p_choice%ssweeps, info)
call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%pvariant, info)
if (psb_toupper(p_choice%solve)=='MUMPS') &
& call prec%set('mumps_loc_glob','local_solver',info)
call prec%set('sub_fillin', p_choice%fill, info)
@ -332,6 +342,8 @@ program amg_d_pde2d
call prec%set('smoother_type', p_choice%smther, info)
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%pvariant, info)
select case (psb_toupper(p_choice%smther))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
@ -362,6 +374,8 @@ program amg_d_pde2d
if (psb_toupper(p_choice%smther2) /= 'NONE') then
call prec%set('smoother_type', p_choice%smther2, info,pos='post')
call prec%set('smoother_sweeps', p_choice%jsweeps2, info,pos='post')
call prec%set('poly_degree', p_choice%degree2, info,pos='post')
call prec%set('poly_variant', p_choice%pvariant2, info,pos='post')
select case (psb_toupper(p_choice%smther2))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
! do nothing
@ -577,6 +591,8 @@ contains
! First smoother / 1-lev preconditioner
call read_data(prec%smther,inp_unit) ! smoother type
call read_data(prec%jsweeps,inp_unit) ! (pre-)smoother / 1-lev prec sweeps
call read_data(prec%degree,inp_unit) ! (pre-)smoother / 1-lev prec sweeps
call read_data(prec%pvariant,inp_unit) !
call read_data(prec%novr,inp_unit) ! number of overlap layers
call read_data(prec%restr,inp_unit) ! restriction over application of AS
call read_data(prec%prol,inp_unit) ! prolongation over application of AS
@ -589,6 +605,8 @@ contains
! Second smoother/ AMG post-smoother (if NONE ignored in main)
call read_data(prec%smther2,inp_unit) ! smoother type
call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps
call read_data(prec%degree2,inp_unit) ! (post-)smoother sweeps
call read_data(prec%pvariant2,inp_unit) !
call read_data(prec%novr2,inp_unit) ! number of overlap layers
call read_data(prec%restr2,inp_unit) ! restriction over application of AS
call read_data(prec%prol2,inp_unit) ! prolongation over application of AS
@ -659,6 +677,8 @@ contains
! broadcast first (pre-)smoother / 1-lev prec data
call psb_bcast(ctxt,prec%smther)
call psb_bcast(ctxt,prec%jsweeps)
call psb_bcast(ctxt,prec%degree)
call psb_bcast(ctxt,prec%pvariant)
call psb_bcast(ctxt,prec%novr)
call psb_bcast(ctxt,prec%restr)
call psb_bcast(ctxt,prec%prol)
@ -671,6 +691,8 @@ contains
! broadcast second (post-)smoother
call psb_bcast(ctxt,prec%smther2)
call psb_bcast(ctxt,prec%jsweeps2)
call psb_bcast(ctxt,prec%degree2)
call psb_bcast(ctxt,prec%pvariant2)
call psb_bcast(ctxt,prec%novr2)
call psb_bcast(ctxt,prec%restr2)
call psb_bcast(ctxt,prec%prol2)

@ -38,52 +38,52 @@ module amg_d_pde2d_base_mod
use psb_base_mod, only : psb_dpk_, dzero, done
real(psb_dpk_), save, private :: epsilon=done/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm2d_base(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm2d_base
!
! functions parametrizing the differential equation
!
function b1(x,y)
function b1_base(x,y)
implicit none
real(psb_dpk_) :: b1
real(psb_dpk_) :: b1_base
real(psb_dpk_), intent(in) :: x,y
b1 = dzero/1.414_psb_dpk_
end function b1
function b2(x,y)
b1_base = dzero/1.414_psb_dpk_
end function b1_base
function b2_base(x,y)
implicit none
real(psb_dpk_) :: b2
real(psb_dpk_) :: b2_base
real(psb_dpk_), intent(in) :: x,y
b2 = dzero/1.414_psb_dpk_
end function b2
function c(x,y)
b2_base = dzero/1.414_psb_dpk_
end function b2_base
function c_base(x,y)
implicit none
real(psb_dpk_) :: c
real(psb_dpk_) :: c_base
real(psb_dpk_), intent(in) :: x,y
c = dzero
end function c
function a1(x,y)
c_base = dzero
end function c_base
function a1_base(x,y)
implicit none
real(psb_dpk_) :: a1
real(psb_dpk_) :: a1_base
real(psb_dpk_), intent(in) :: x,y
a1=done*epsilon
end function a1
function a2(x,y)
a1_base=done*epsilon
end function a1_base
function a2_base(x,y)
implicit none
real(psb_dpk_) :: a2
real(psb_dpk_) :: a2_base
real(psb_dpk_), intent(in) :: x,y
a2=done*epsilon
end function a2
function g(x,y)
a2_base=done*epsilon
end function a2_base
function g_base(x,y)
implicit none
real(psb_dpk_) :: g
real(psb_dpk_) :: g_base
real(psb_dpk_), intent(in) :: x,y
g = dzero
g_base = dzero
if (x == done) then
g = done
g_base = done
else if (x == dzero) then
g = done
g_base = done
end if
end function g
end function g_base
end module amg_d_pde2d_base_mod

@ -38,10 +38,10 @@ module amg_d_pde2d_box_mod
use psb_base_mod, only : psb_dpk_, dzero, done
real(psb_dpk_), save, private :: epsilon=done/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm2d_box(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm2d_box
!
! functions parametrizing the differential equation
!

@ -38,10 +38,10 @@ module amg_d_pde2d_exp_mod
use psb_base_mod, only : psb_dpk_, done, dzero
real(psb_dpk_), save, private :: epsilon=done/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm2d_exp(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm2d_exp
!
! functions parametrizing the differential equation
!

@ -0,0 +1,89 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module amg_d_pde2d_gauss_mod
use psb_base_mod, only : psb_dpk_, done, dzero
real(psb_dpk_), save, private :: epsilon=done/80
contains
subroutine pde_set_parm2d_gauss(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm2d_gauss
!
! functions parametrizing the differential equation
!
function b1_gauss(x,y)
implicit none
real(psb_dpk_) :: b1_gauss
real(psb_dpk_), intent(in) :: x,y
b1_gauss=done/sqrt(3.0_psb_dpk_)-2*x*exp(-(x**2+y**2))
end function b1_gauss
function b2_gauss(x,y)
implicit none
real(psb_dpk_) :: b2_gauss
real(psb_dpk_), intent(in) :: x,y
b2_gauss=done/sqrt(3.0_psb_dpk_)-2*y*exp(-(x**2+y**2))
end function b2_gauss
function c_gauss(x,y)
implicit none
real(psb_dpk_) :: c_gauss
real(psb_dpk_), intent(in) :: x,y
c_gauss=dzero
end function c_gauss
function a1_gauss(x,y)
implicit none
real(psb_dpk_) :: a1_gauss
real(psb_dpk_), intent(in) :: x,y
a1_gauss=epsilon*exp(-(x**2+y**2))
end function a1_gauss
function a2_gauss(x,y)
implicit none
real(psb_dpk_) :: a2_gauss
real(psb_dpk_), intent(in) :: x,y
a2_gauss=epsilon*exp(-(x**2+y**2))
end function a2_gauss
function g_gauss(x,y)
implicit none
real(psb_dpk_) :: g_gauss
real(psb_dpk_), intent(in) :: x,y
g_gauss = dzero
if (x == done) then
g_gauss = done
else if (x == dzero) then
g_gauss = done
end if
end function g_gauss
end module amg_d_pde2d_gauss_mod

@ -72,6 +72,7 @@ program amg_d_pde3d
use data_input
use amg_d_pde3d_base_mod
use amg_d_pde3d_exp_mod
use amg_d_pde3d_box_mod
use amg_d_pde3d_gauss_mod
use amg_d_genpde_mod
#if defined(OPENMP)
@ -201,7 +202,7 @@ program amg_d_pde3d
! other variables
integer(psb_ipk_) :: info, i, k
character(len=20) :: name,ch_err
type(psb_d_csr_sparse_mat) :: amold
info=psb_success_
@ -247,10 +248,13 @@ program amg_d_pde3d
select case(psb_toupper(trim(pdecoeff)))
case("CONST")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info)
& a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info)
case("EXP")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info)
case("BOX")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_box,a2_box,a3_box,b1_box,b2_box,b3_box,c_box,g_box,info)
case("GAUSS")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
@ -412,7 +416,6 @@ program amg_d_pde3d
call prec%set('coarse_sweeps', p_choice%cjswp, info)
end select
!!$ call prec%descr(info,iout=psb_out_unit)
! build the preconditioner
call psb_barrier(ctxt)
@ -425,7 +428,7 @@ program amg_d_pde3d
end if
call psb_barrier(ctxt)
t1 = psb_wtime()
call prec%smoothers_build(a,desc_a,info,amold=amold)
call prec%smoothers_build(a,desc_a,info)
tprec = psb_wtime()-t1
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_smoothers_bld')

@ -38,64 +38,64 @@ module amg_d_pde3d_base_mod
use psb_base_mod, only : psb_dpk_, done, dzero
real(psb_dpk_), save, private :: epsilon=done/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm3d_base(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm3d_base
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
function b1_base(x,y,z)
implicit none
real(psb_dpk_) :: b1
real(psb_dpk_) :: b1_base
real(psb_dpk_), intent(in) :: x,y,z
b1=dzero/sqrt(3.0_psb_dpk_)
end function b1
function b2(x,y,z)
b1_base=dzero/sqrt(3.0_psb_dpk_)
end function b1_base
function b2_base(x,y,z)
implicit none
real(psb_dpk_) :: b2
real(psb_dpk_) :: b2_base
real(psb_dpk_), intent(in) :: x,y,z
b2=dzero/sqrt(3.0_psb_dpk_)
end function b2
function b3(x,y,z)
b2_base=dzero/sqrt(3.0_psb_dpk_)
end function b2_base
function b3_base(x,y,z)
implicit none
real(psb_dpk_) :: b3
real(psb_dpk_) :: b3_base
real(psb_dpk_), intent(in) :: x,y,z
b3=dzero/sqrt(3.0_psb_dpk_)
end function b3
function c(x,y,z)
b3_base=dzero/sqrt(3.0_psb_dpk_)
end function b3_base
function c_base(x,y,z)
implicit none
real(psb_dpk_) :: c
real(psb_dpk_) :: c_base
real(psb_dpk_), intent(in) :: x,y,z
c=dzero
end function c
function a1(x,y,z)
c_base=dzero
end function c_base
function a1_base(x,y,z)
implicit none
real(psb_dpk_) :: a1
real(psb_dpk_) :: a1_base
real(psb_dpk_), intent(in) :: x,y,z
a1=epsilon
end function a1
function a2(x,y,z)
a1_base=epsilon
end function a1_base
function a2_base(x,y,z)
implicit none
real(psb_dpk_) :: a2
real(psb_dpk_) :: a2_base
real(psb_dpk_), intent(in) :: x,y,z
a2=epsilon
end function a2
function a3(x,y,z)
a2_base=epsilon
end function a2_base
function a3_base(x,y,z)
implicit none
real(psb_dpk_) :: a3
real(psb_dpk_) :: a3_base
real(psb_dpk_), intent(in) :: x,y,z
a3=epsilon
end function a3
function g(x,y,z)
a3_base=epsilon
end function a3_base
function g_base(x,y,z)
implicit none
real(psb_dpk_) :: g
real(psb_dpk_) :: g_base
real(psb_dpk_), intent(in) :: x,y,z
g = dzero
g_base = dzero
if (x == done) then
g = done
g_base = done
else if (x == dzero) then
g = done
g_base = done
end if
end function g
end function g_base
end module amg_d_pde3d_base_mod

@ -0,0 +1,101 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module amg_d_pde3d_box_mod
use psb_base_mod, only : psb_dpk_, done, dzero
real(psb_dpk_), save, private :: epsilon=done/80
contains
subroutine pde_set_parm3d_box(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm3d_box
!
! functions parametrizing the differential equation
!
function b1_box(x,y,z)
implicit none
real(psb_dpk_) :: b1_box
real(psb_dpk_), intent(in) :: x,y,z
b1_box=done/sqrt(3.0_psb_dpk_)
end function b1_box
function b2_box(x,y,z)
implicit none
real(psb_dpk_) :: b2_box
real(psb_dpk_), intent(in) :: x,y,z
b2_box=done/sqrt(3.0_psb_dpk_)
end function b2_box
function b3_box(x,y,z)
implicit none
real(psb_dpk_) :: b3_box
real(psb_dpk_), intent(in) :: x,y,z
b3_box=done/sqrt(3.0_psb_dpk_)
end function b3_box
function c_box(x,y,z)
implicit none
real(psb_dpk_) :: c_box
real(psb_dpk_), intent(in) :: x,y,z
c_box=dzero
end function c_box
function a1_box(x,y,z)
implicit none
real(psb_dpk_) :: a1_box
real(psb_dpk_), intent(in) :: x,y,z
a1_box=epsilon
end function a1_box
function a2_box(x,y,z)
implicit none
real(psb_dpk_) :: a2_box
real(psb_dpk_), intent(in) :: x,y,z
a2_box=epsilon
end function a2_box
function a3_box(x,y,z)
implicit none
real(psb_dpk_) :: a3_box
real(psb_dpk_), intent(in) :: x,y,z
a3_box=epsilon
end function a3_box
function g_box(x,y,z)
implicit none
real(psb_dpk_) :: g_box
real(psb_dpk_), intent(in) :: x,y,z
g_box= dzero
if (x == done) then
g_box = done
else if (x == dzero) then
g_box = done
end if
end function g_box
end module amg_d_pde3d_box_mod

@ -38,10 +38,10 @@ module amg_d_pde3d_exp_mod
use psb_base_mod, only : psb_dpk_, done, dzero
real(psb_dpk_), save, private :: epsilon=done/160
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm3d_exp(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm3d_exp
!
! functions parametrizing the differential equation
!

@ -38,10 +38,10 @@ module amg_d_pde3d_gauss_mod
use psb_base_mod, only : psb_dpk_, done, dzero
real(psb_dpk_), save, private :: epsilon=done/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm3d_gauss(dat)
real(psb_dpk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm3d_gauss
!
! functions parametrizing the differential equation
!

@ -72,6 +72,7 @@ program amg_s_pde2d
use amg_s_pde2d_base_mod
use amg_s_pde2d_exp_mod
use amg_s_pde2d_box_mod
use amg_s_pde2d_gauss_mod
use amg_s_genpde_mod
#if defined(OPENMP)
use omp_lib
@ -124,16 +125,16 @@ program amg_s_pde2d
integer(psb_ipk_) :: outer_sweeps ! number of outer sweeps: sweeps for 1-level,
! AMG cycles for ML
! general AMG data
character(len=16) :: mlcycle ! AMG cycle type
character(len=32) :: mlcycle ! AMG cycle type
integer(psb_ipk_) :: maxlevs ! maximum number of levels in AMG preconditioner
! AMG aggregation
character(len=16) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED
character(len=16) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC
character(len=16) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP
character(len=32) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED
character(len=32) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC
character(len=32) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP
integer(psb_ipk_) :: aggr_size ! Requested size of the aggregates for MATCHBOXP
character(len=16) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE
character(len=16) :: aggr_filter ! filtering: FILTER, NO_FILTER
character(len=32) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE
character(len=32) :: aggr_filter ! filtering: FILTER, NO_FILTER
real(psb_spk_) :: mncrratio ! minimum aggregation ratio
real(psb_spk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector
integer(psb_ipk_) :: thrvsz ! size of threshold vector
@ -141,39 +142,43 @@ program amg_s_pde2d
integer(psb_ipk_) :: csizepp ! minimum size of coarsest matrix per process
! AMG smoother or pre-smoother; also 1-lev preconditioner
character(len=16) :: smther ! (pre-)smoother type: BJAC, AS
character(len=32) :: smther ! (pre-)smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps
integer(psb_ipk_) :: degree ! degree for polynomial smoother
character(len=32) :: pvariant ! polynomial variant
integer(psb_ipk_) :: novr ! number of overlap layers
character(len=16) :: restr ! restriction over application of AS
character(len=16) :: prol ! prolongation over application of AS
character(len=16) :: solve ! local subsolver type: ILU, MILU, ILUT,
character(len=32) :: restr ! restriction over application of AS
character(len=32) :: prol ! prolongation over application of AS
character(len=32) :: solve ! local subsolver type: ILU, MILU, ILUT,
! UMF, MUMPS, SLU, FWGS, BWGS, JAC
integer(psb_ipk_) :: ssweeps ! inner solver sweeps
character(len=16) :: variant ! AINV variant: LLK, etc
character(len=32) :: variant ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill ! Inverse fill-in for INVK
real(psb_spk_) :: thr ! threshold for ILUT factorization
! AMG post-smoother; ignored by 1-lev preconditioner
character(len=16) :: smther2 ! post-smoother type: BJAC, AS
character(len=32) :: smther2 ! post-smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps
integer(psb_ipk_) :: degree2 ! degree for polynomial smoother
character(len=32) :: pvariant2 ! polynomial variant
integer(psb_ipk_) :: novr2 ! number of overlap layers
character(len=16) :: restr2 ! restriction over application of AS
character(len=16) :: prol2 ! prolongation over application of AS
character(len=16) :: solve2 ! local subsolver type: ILU, MILU, ILUT,
character(len=32) :: restr2 ! restriction over application of AS
character(len=32) :: prol2 ! prolongation over application of AS
character(len=32) :: solve2 ! local subsolver type: ILU, MILU, ILUT,
! UMF, MUMPS, SLU, FWGS, BWGS, JAC
integer(psb_ipk_) :: ssweeps2 ! inner solver sweeps
character(len=16) :: variant2 ! AINV variant: LLK, etc
character(len=32) :: variant2 ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill2 ! Inverse fill-in for INVK
real(psb_spk_) :: thr2 ! threshold for ILUT factorization
! coarsest-level solver
character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST
character(len=16) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr.
character(len=32) :: cmat ! coarsest matrix layout: REPL, DIST
character(len=32) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr.
! mat.); UMF, MUMPS, SLU, ILU, ILUT, MILU
! (repl. mat.)
character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT,
character(len=32) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT,
! MILU, UMF, MUMPS, SLU
integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization
real(psb_spk_) :: cthres ! threshold for ILUT factorization
@ -241,13 +246,16 @@ program amg_s_pde2d
select case(psb_toupper(trim(pdecoeff)))
case("CONST")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,b1,b2,c,g,info)
& a1_base,a2_base,b1_base,b2_base,c_base,g_base,info)
case("EXP")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_exp,a2_exp,b1_exp,b2_exp,c_exp,g_exp,info)
case("BOX")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_box,a2_box,b1_box,b2_box,c_box,g_box,info)
case("GAUSS")
call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_gauss,a2_gauss,b1_gauss,b2_gauss,c_gauss,g_gauss,info)
case default
info=psb_err_from_subroutine_
ch_err='amg_gen_pdecoeff'
@ -281,10 +289,12 @@ program amg_s_pde2d
! 1-level sweeps from "outer_sweeps"
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
case ('BJAC')
case ('BJAC','POLY')
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('sub_solve', p_choice%solve, info)
call prec%set('solver_sweeps', p_choice%ssweeps, info)
call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%pvariant, info)
if (psb_toupper(p_choice%solve)=='MUMPS') &
& call prec%set('mumps_loc_glob','local_solver',info)
call prec%set('sub_fillin', p_choice%fill, info)
@ -332,6 +342,8 @@ program amg_s_pde2d
call prec%set('smoother_type', p_choice%smther, info)
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%pvariant, info)
select case (psb_toupper(p_choice%smther))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
@ -362,6 +374,8 @@ program amg_s_pde2d
if (psb_toupper(p_choice%smther2) /= 'NONE') then
call prec%set('smoother_type', p_choice%smther2, info,pos='post')
call prec%set('smoother_sweeps', p_choice%jsweeps2, info,pos='post')
call prec%set('poly_degree', p_choice%degree2, info,pos='post')
call prec%set('poly_variant', p_choice%pvariant2, info,pos='post')
select case (psb_toupper(p_choice%smther2))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
! do nothing
@ -577,6 +591,8 @@ contains
! First smoother / 1-lev preconditioner
call read_data(prec%smther,inp_unit) ! smoother type
call read_data(prec%jsweeps,inp_unit) ! (pre-)smoother / 1-lev prec sweeps
call read_data(prec%degree,inp_unit) ! (pre-)smoother / 1-lev prec sweeps
call read_data(prec%pvariant,inp_unit) !
call read_data(prec%novr,inp_unit) ! number of overlap layers
call read_data(prec%restr,inp_unit) ! restriction over application of AS
call read_data(prec%prol,inp_unit) ! prolongation over application of AS
@ -589,6 +605,8 @@ contains
! Second smoother/ AMG post-smoother (if NONE ignored in main)
call read_data(prec%smther2,inp_unit) ! smoother type
call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps
call read_data(prec%degree2,inp_unit) ! (post-)smoother sweeps
call read_data(prec%pvariant2,inp_unit) !
call read_data(prec%novr2,inp_unit) ! number of overlap layers
call read_data(prec%restr2,inp_unit) ! restriction over application of AS
call read_data(prec%prol2,inp_unit) ! prolongation over application of AS
@ -659,6 +677,8 @@ contains
! broadcast first (pre-)smoother / 1-lev prec data
call psb_bcast(ctxt,prec%smther)
call psb_bcast(ctxt,prec%jsweeps)
call psb_bcast(ctxt,prec%degree)
call psb_bcast(ctxt,prec%pvariant)
call psb_bcast(ctxt,prec%novr)
call psb_bcast(ctxt,prec%restr)
call psb_bcast(ctxt,prec%prol)
@ -671,6 +691,8 @@ contains
! broadcast second (post-)smoother
call psb_bcast(ctxt,prec%smther2)
call psb_bcast(ctxt,prec%jsweeps2)
call psb_bcast(ctxt,prec%degree2)
call psb_bcast(ctxt,prec%pvariant2)
call psb_bcast(ctxt,prec%novr2)
call psb_bcast(ctxt,prec%restr2)
call psb_bcast(ctxt,prec%prol2)

@ -38,52 +38,52 @@ module amg_s_pde2d_base_mod
use psb_base_mod, only : psb_spk_, szero, sone
real(psb_spk_), save, private :: epsilon=sone/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm2d_base(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm2d_base
!
! functions parametrizing the differential equation
!
function b1(x,y)
function b1_base(x,y)
implicit none
real(psb_spk_) :: b1
real(psb_spk_) :: b1_base
real(psb_spk_), intent(in) :: x,y
b1 = szero/1.414_psb_spk_
end function b1
function b2(x,y)
b1_base = szero/1.414_psb_spk_
end function b1_base
function b2_base(x,y)
implicit none
real(psb_spk_) :: b2
real(psb_spk_) :: b2_base
real(psb_spk_), intent(in) :: x,y
b2 = szero/1.414_psb_spk_
end function b2
function c(x,y)
b2_base = szero/1.414_psb_spk_
end function b2_base
function c_base(x,y)
implicit none
real(psb_spk_) :: c
real(psb_spk_) :: c_base
real(psb_spk_), intent(in) :: x,y
c = szero
end function c
function a1(x,y)
c_base = szero
end function c_base
function a1_base(x,y)
implicit none
real(psb_spk_) :: a1
real(psb_spk_) :: a1_base
real(psb_spk_), intent(in) :: x,y
a1=sone*epsilon
end function a1
function a2(x,y)
a1_base=sone*epsilon
end function a1_base
function a2_base(x,y)
implicit none
real(psb_spk_) :: a2
real(psb_spk_) :: a2_base
real(psb_spk_), intent(in) :: x,y
a2=sone*epsilon
end function a2
function g(x,y)
a2_base=sone*epsilon
end function a2_base
function g_base(x,y)
implicit none
real(psb_spk_) :: g
real(psb_spk_) :: g_base
real(psb_spk_), intent(in) :: x,y
g = szero
g_base = szero
if (x == sone) then
g = sone
g_base = sone
else if (x == szero) then
g = sone
g_base = sone
end if
end function g
end function g_base
end module amg_s_pde2d_base_mod

@ -38,10 +38,10 @@ module amg_s_pde2d_box_mod
use psb_base_mod, only : psb_spk_, szero, sone
real(psb_spk_), save, private :: epsilon=sone/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm2d_box(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm2d_box
!
! functions parametrizing the differential equation
!

@ -38,10 +38,10 @@ module amg_s_pde2d_exp_mod
use psb_base_mod, only : psb_spk_, sone, szero
real(psb_spk_), save, private :: epsilon=sone/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm2d_exp(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm2d_exp
!
! functions parametrizing the differential equation
!

@ -0,0 +1,89 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module amg_s_pde2d_gauss_mod
use psb_base_mod, only : psb_spk_, sone, szero
real(psb_spk_), save, private :: epsilon=sone/80
contains
subroutine pde_set_parm2d_gauss(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm2d_gauss
!
! functions parametrizing the differential equation
!
function b1_gauss(x,y)
implicit none
real(psb_spk_) :: b1_gauss
real(psb_spk_), intent(in) :: x,y
b1_gauss=sone/sqrt(3.0_psb_spk_)-2*x*exp(-(x**2+y**2))
end function b1_gauss
function b2_gauss(x,y)
implicit none
real(psb_spk_) :: b2_gauss
real(psb_spk_), intent(in) :: x,y
b2_gauss=sone/sqrt(3.0_psb_spk_)-2*y*exp(-(x**2+y**2))
end function b2_gauss
function c_gauss(x,y)
implicit none
real(psb_spk_) :: c_gauss
real(psb_spk_), intent(in) :: x,y
c_gauss=szero
end function c_gauss
function a1_gauss(x,y)
implicit none
real(psb_spk_) :: a1_gauss
real(psb_spk_), intent(in) :: x,y
a1_gauss=epsilon*exp(-(x**2+y**2))
end function a1_gauss
function a2_gauss(x,y)
implicit none
real(psb_spk_) :: a2_gauss
real(psb_spk_), intent(in) :: x,y
a2_gauss=epsilon*exp(-(x**2+y**2))
end function a2_gauss
function g_gauss(x,y)
implicit none
real(psb_spk_) :: g_gauss
real(psb_spk_), intent(in) :: x,y
g_gauss = szero
if (x == sone) then
g_gauss = sone
else if (x == szero) then
g_gauss = sone
end if
end function g_gauss
end module amg_s_pde2d_gauss_mod

@ -72,6 +72,7 @@ program amg_s_pde3d
use data_input
use amg_s_pde3d_base_mod
use amg_s_pde3d_exp_mod
use amg_s_pde3d_box_mod
use amg_s_pde3d_gauss_mod
use amg_s_genpde_mod
#if defined(OPENMP)
@ -125,16 +126,16 @@ program amg_s_pde3d
integer(psb_ipk_) :: outer_sweeps ! number of outer sweeps: sweeps for 1-level,
! AMG cycles for ML
! general AMG data
character(len=16) :: mlcycle ! AMG cycle type
character(len=32) :: mlcycle ! AMG cycle type
integer(psb_ipk_) :: maxlevs ! maximum number of levels in AMG preconditioner
! AMG aggregation
character(len=16) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED
character(len=16) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC
character(len=16) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP
character(len=32) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED
character(len=32) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC
character(len=32) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP
integer(psb_ipk_) :: aggr_size ! Requested size of the aggregates for MATCHBOXP
character(len=16) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE
character(len=16) :: aggr_filter ! filtering: FILTER, NO_FILTER
character(len=32) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE
character(len=32) :: aggr_filter ! filtering: FILTER, NO_FILTER
real(psb_spk_) :: mncrratio ! minimum aggregation ratio
real(psb_spk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector
integer(psb_ipk_) :: thrvsz ! size of threshold vector
@ -142,39 +143,43 @@ program amg_s_pde3d
integer(psb_ipk_) :: csizepp ! minimum size of coarsest matrix per process
! AMG smoother or pre-smoother; also 1-lev preconditioner
character(len=16) :: smther ! (pre-)smoother type: BJAC, AS
character(len=32) :: smther ! (pre-)smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps
integer(psb_ipk_) :: degree ! degree for polynomial smoother
character(len=32) :: pvariant ! polynomial variant
integer(psb_ipk_) :: novr ! number of overlap layers
character(len=16) :: restr ! restriction over application of AS
character(len=16) :: prol ! prolongation over application of AS
character(len=16) :: solve ! local subsolver type: ILU, MILU, ILUT,
character(len=32) :: restr ! restriction over application of AS
character(len=32) :: prol ! prolongation over application of AS
character(len=32) :: solve ! local subsolver type: ILU, MILU, ILUT,
! UMF, MUMPS, SLU, FWGS, BWGS, JAC
integer(psb_ipk_) :: ssweeps ! inner solver sweeps
character(len=16) :: variant ! AINV variant: LLK, etc
character(len=32) :: variant ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill ! Inverse fill-in for INVK
real(psb_spk_) :: thr ! threshold for ILUT factorization
! AMG post-smoother; ignored by 1-lev preconditioner
character(len=16) :: smther2 ! post-smoother type: BJAC, AS
character(len=32) :: smther2 ! post-smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps
integer(psb_ipk_) :: degree2 ! degree for polynomial smoother
character(len=32) :: pvariant2 ! polynomial variant
integer(psb_ipk_) :: novr2 ! number of overlap layers
character(len=16) :: restr2 ! restriction over application of AS
character(len=16) :: prol2 ! prolongation over application of AS
character(len=16) :: solve2 ! local subsolver type: ILU, MILU, ILUT,
character(len=32) :: restr2 ! restriction over application of AS
character(len=32) :: prol2 ! prolongation over application of AS
character(len=32) :: solve2 ! local subsolver type: ILU, MILU, ILUT,
! UMF, MUMPS, SLU, FWGS, BWGS, JAC
integer(psb_ipk_) :: ssweeps2 ! inner solver sweeps
character(len=16) :: variant2 ! AINV variant: LLK, etc
character(len=32) :: variant2 ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill2 ! Inverse fill-in for INVK
real(psb_spk_) :: thr2 ! threshold for ILUT factorization
! coarsest-level solver
character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST
character(len=16) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr.
character(len=32) :: cmat ! coarsest matrix layout: REPL, DIST
character(len=32) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr.
! mat.); UMF, MUMPS, SLU, ILU, ILUT, MILU
! (repl. mat.)
character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT,
character(len=32) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT,
! MILU, UMF, MUMPS, SLU
integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization
real(psb_spk_) :: cthres ! threshold for ILUT factorization
@ -197,7 +202,7 @@ program amg_s_pde3d
! other variables
integer(psb_ipk_) :: info, i, k
character(len=20) :: name,ch_err
type(psb_s_csr_sparse_mat) :: amold
info=psb_success_
@ -243,10 +248,13 @@ program amg_s_pde3d
select case(psb_toupper(trim(pdecoeff)))
case("CONST")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info)
& a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info)
case("EXP")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info)
case("BOX")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_box,a2_box,a3_box,b1_box,b2_box,b3_box,c_box,g_box,info)
case("GAUSS")
call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
& a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
@ -285,10 +293,12 @@ program amg_s_pde3d
! 1-level sweeps from "outer_sweeps"
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
case ('BJAC')
case ('BJAC','POLY')
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('sub_solve', p_choice%solve, info)
call prec%set('solver_sweeps', p_choice%ssweeps, info)
call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%pvariant, info)
if (psb_toupper(p_choice%solve)=='MUMPS') &
& call prec%set('mumps_loc_glob','local_solver',info)
call prec%set('sub_fillin', p_choice%fill, info)
@ -336,6 +346,8 @@ program amg_s_pde3d
call prec%set('smoother_type', p_choice%smther, info)
call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%pvariant, info)
select case (psb_toupper(p_choice%smther))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
@ -366,6 +378,8 @@ program amg_s_pde3d
if (psb_toupper(p_choice%smther2) /= 'NONE') then
call prec%set('smoother_type', p_choice%smther2, info,pos='post')
call prec%set('smoother_sweeps', p_choice%jsweeps2, info,pos='post')
call prec%set('poly_degree', p_choice%degree2, info,pos='post')
call prec%set('poly_variant', p_choice%pvariant2, info,pos='post')
select case (psb_toupper(p_choice%smther2))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
! do nothing
@ -414,7 +428,7 @@ program amg_s_pde3d
end if
call psb_barrier(ctxt)
t1 = psb_wtime()
call prec%smoothers_build(a,desc_a,info,amold=amold)
call prec%smoothers_build(a,desc_a,info)
tprec = psb_wtime()-t1
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_smoothers_bld')
@ -581,6 +595,8 @@ contains
! First smoother / 1-lev preconditioner
call read_data(prec%smther,inp_unit) ! smoother type
call read_data(prec%jsweeps,inp_unit) ! (pre-)smoother / 1-lev prec sweeps
call read_data(prec%degree,inp_unit) ! (pre-)smoother / 1-lev prec sweeps
call read_data(prec%pvariant,inp_unit) !
call read_data(prec%novr,inp_unit) ! number of overlap layers
call read_data(prec%restr,inp_unit) ! restriction over application of AS
call read_data(prec%prol,inp_unit) ! prolongation over application of AS
@ -593,6 +609,8 @@ contains
! Second smoother/ AMG post-smoother (if NONE ignored in main)
call read_data(prec%smther2,inp_unit) ! smoother type
call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps
call read_data(prec%degree2,inp_unit) ! (post-)smoother sweeps
call read_data(prec%pvariant2,inp_unit) !
call read_data(prec%novr2,inp_unit) ! number of overlap layers
call read_data(prec%restr2,inp_unit) ! restriction over application of AS
call read_data(prec%prol2,inp_unit) ! prolongation over application of AS
@ -663,6 +681,8 @@ contains
! broadcast first (pre-)smoother / 1-lev prec data
call psb_bcast(ctxt,prec%smther)
call psb_bcast(ctxt,prec%jsweeps)
call psb_bcast(ctxt,prec%degree)
call psb_bcast(ctxt,prec%pvariant)
call psb_bcast(ctxt,prec%novr)
call psb_bcast(ctxt,prec%restr)
call psb_bcast(ctxt,prec%prol)
@ -675,6 +695,8 @@ contains
! broadcast second (post-)smoother
call psb_bcast(ctxt,prec%smther2)
call psb_bcast(ctxt,prec%jsweeps2)
call psb_bcast(ctxt,prec%degree2)
call psb_bcast(ctxt,prec%pvariant2)
call psb_bcast(ctxt,prec%novr2)
call psb_bcast(ctxt,prec%restr2)
call psb_bcast(ctxt,prec%prol2)

@ -38,64 +38,64 @@ module amg_s_pde3d_base_mod
use psb_base_mod, only : psb_spk_, sone, szero
real(psb_spk_), save, private :: epsilon=sone/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm3d_base(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm3d_base
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
function b1_base(x,y,z)
implicit none
real(psb_spk_) :: b1
real(psb_spk_) :: b1_base
real(psb_spk_), intent(in) :: x,y,z
b1=sone/sqrt(3.0_psb_spk_)
end function b1
function b2(x,y,z)
b1_base=szero/sqrt(3.0_psb_spk_)
end function b1_base
function b2_base(x,y,z)
implicit none
real(psb_spk_) :: b2
real(psb_spk_) :: b2_base
real(psb_spk_), intent(in) :: x,y,z
b2=sone/sqrt(3.0_psb_spk_)
end function b2
function b3(x,y,z)
b2_base=szero/sqrt(3.0_psb_spk_)
end function b2_base
function b3_base(x,y,z)
implicit none
real(psb_spk_) :: b3
real(psb_spk_) :: b3_base
real(psb_spk_), intent(in) :: x,y,z
b3=sone/sqrt(3.0_psb_spk_)
end function b3
function c(x,y,z)
b3_base=szero/sqrt(3.0_psb_spk_)
end function b3_base
function c_base(x,y,z)
implicit none
real(psb_spk_) :: c
real(psb_spk_) :: c_base
real(psb_spk_), intent(in) :: x,y,z
c=szero
end function c
function a1(x,y,z)
c_base=szero
end function c_base
function a1_base(x,y,z)
implicit none
real(psb_spk_) :: a1
real(psb_spk_) :: a1_base
real(psb_spk_), intent(in) :: x,y,z
a1=epsilon
end function a1
function a2(x,y,z)
a1_base=epsilon
end function a1_base
function a2_base(x,y,z)
implicit none
real(psb_spk_) :: a2
real(psb_spk_) :: a2_base
real(psb_spk_), intent(in) :: x,y,z
a2=epsilon
end function a2
function a3(x,y,z)
a2_base=epsilon
end function a2_base
function a3_base(x,y,z)
implicit none
real(psb_spk_) :: a3
real(psb_spk_) :: a3_base
real(psb_spk_), intent(in) :: x,y,z
a3=epsilon
end function a3
function g(x,y,z)
a3_base=epsilon
end function a3_base
function g_base(x,y,z)
implicit none
real(psb_spk_) :: g
real(psb_spk_) :: g_base
real(psb_spk_), intent(in) :: x,y,z
g = szero
g_base = szero
if (x == sone) then
g = sone
g_base = sone
else if (x == szero) then
g = sone
g_base = sone
end if
end function g
end function g_base
end module amg_s_pde3d_base_mod

@ -0,0 +1,101 @@
!
!
! AMG4PSBLAS version 1.0
! Algebraic Multigrid Package
! based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
! (C) Copyright 2021
!
! Salvatore Filippone
! Pasqua D'Ambra
! Fabio Durastante
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the AMG4PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module amg_s_pde3d_box_mod
use psb_base_mod, only : psb_spk_, sone, szero
real(psb_spk_), save, private :: epsilon=sone/80
contains
subroutine pde_set_parm3d_box(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm3d_box
!
! functions parametrizing the differential equation
!
function b1_box(x,y,z)
implicit none
real(psb_spk_) :: b1_box
real(psb_spk_), intent(in) :: x,y,z
b1_box=sone/sqrt(3.0_psb_spk_)
end function b1_box
function b2_box(x,y,z)
implicit none
real(psb_spk_) :: b2_box
real(psb_spk_), intent(in) :: x,y,z
b2_box=sone/sqrt(3.0_psb_spk_)
end function b2_box
function b3_box(x,y,z)
implicit none
real(psb_spk_) :: b3_box
real(psb_spk_), intent(in) :: x,y,z
b3_box=sone/sqrt(3.0_psb_spk_)
end function b3_box
function c_box(x,y,z)
implicit none
real(psb_spk_) :: c_box
real(psb_spk_), intent(in) :: x,y,z
c_box=szero
end function c_box
function a1_box(x,y,z)
implicit none
real(psb_spk_) :: a1_box
real(psb_spk_), intent(in) :: x,y,z
a1_box=epsilon
end function a1_box
function a2_box(x,y,z)
implicit none
real(psb_spk_) :: a2_box
real(psb_spk_), intent(in) :: x,y,z
a2_box=epsilon
end function a2_box
function a3_box(x,y,z)
implicit none
real(psb_spk_) :: a3_box
real(psb_spk_), intent(in) :: x,y,z
a3_box=epsilon
end function a3_box
function g_box(x,y,z)
implicit none
real(psb_spk_) :: g_box
real(psb_spk_), intent(in) :: x,y,z
g_box= szero
if (x == sone) then
g_box = sone
else if (x == szero) then
g_box = sone
end if
end function g_box
end module amg_s_pde3d_box_mod

@ -38,10 +38,10 @@ module amg_s_pde3d_exp_mod
use psb_base_mod, only : psb_spk_, sone, szero
real(psb_spk_), save, private :: epsilon=sone/160
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm3d_exp(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm3d_exp
!
! functions parametrizing the differential equation
!

@ -38,10 +38,10 @@ module amg_s_pde3d_gauss_mod
use psb_base_mod, only : psb_spk_, sone, szero
real(psb_spk_), save, private :: epsilon=sone/80
contains
subroutine pde_set_parm(dat)
subroutine pde_set_parm3d_gauss(dat)
real(psb_spk_), intent(in) :: dat
epsilon = dat
end subroutine pde_set_parm
end subroutine pde_set_parm3d_gauss
!
! functions parametrizing the differential equation
!

@ -1,7 +1,7 @@
%%%%%%%%%%% General arguments % Lines starting with % are ignored.
CSR ! Storage format CSR COO JAD
0200 ! IDIM; domain size. Linear system size is IDIM**2
CONST ! PDECOEFF: CONST, EXP, BOX Coefficients of the PDE
CONST ! PDECOEFF: CONST, EXP, BOX, GAUSS Coefficients of the PDE
CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
2 ! ISTOPC
00500 ! ITMAX
@ -9,23 +9,33 @@ CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS F
30 ! IRST (restart for RGMRES and BiCGSTABL)
1.d-6 ! EPS
%%%%%%%%%%% Main preconditioner choices %%%%%%%%%%%%%%%%
ML-VCYCLE-BJAC-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars)
ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML
ML-VBM-VCYCLE-FBGS-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars)
ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY
%
%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
BJAC ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous.
1 ! Number of sweeps for smoother
FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r 1-level, repeats previous.
6 ! Number of sweeps for smoother
1 ! degree for polynomial smoother
POLY_LOTTES_BETA ! Polynomial variant
% Fields to be added for POLY
% POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER
% POLY_RHO_ESTIMATE_ITERATIONS default = 20
% POLY_RHO_BA set to value
%
0 ! Number of overlap layers for AS preconditioner
HALO ! AS restriction operator: NONE HALO
NONE ! AS prolongation operator: NONE SUM AVG
ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF
8 ! Inner solver sweeps (GS and JACOBI)
LLK ! AINV variant, ignored otherwise
L1-JACOBI ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF
1 ! Inner solver sweeps (GS and JACOBI)
LLK ! AINV variant
0 ! Fill level P for ILU(P) and ILU(T,P)
1 ! Inverse Fill level P for INVK
1.d-4 ! Threshold T for ILU(T,P)
%%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%%
NONE ! Second (post) smoother, ignored if NONE
1 ! Number of sweeps for (post) smoother
6 ! Number of sweeps for (post) smoother
1 ! degree for polynomial smoother
POLY_LOTTES_BETA ! Polynomial variant
0 ! Number of overlap layers for AS preconditioner
HALO ! AS restriction operator: NONE HALO
NONE ! AS prolongation operator: NONE SUM AVG

@ -1,7 +1,7 @@
%%%%%%%%%%% General arguments % Lines starting with % are ignored.
CSR ! Storage format CSR COO JAD
0150 ! IDIM; domain size. Linear system size is IDIM**3
CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE
CONST ! PDECOEFF: CONST, EXP, BOX, GAUSS Coefficients of the PDE
CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
2 ! ISTOPC
00500 ! ITMAX
@ -12,16 +12,16 @@ CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS F
ML-VBM-VCYCLE-FBGS-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars)
ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY
%
%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r 1-level, repeats previous.
6 ! Number of sweeps for smoother
1 ! degree for polynomial smoother
POLY_LOTTES_BETA ! Polynomial variant
% Fields to be added for POLY
% POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER
% POLY_RHO_ESTIMATE_ITERATIONS default = 20
% POLY_RHO_BA set to value
%
%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
L1-JACOBI ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r 1-level, repeats previous.
6 ! Number of sweeps for smoother
1 ! degree for polynomial smoother
POLY_LOTTES_BETA ! Polynomial variant
0 ! Number of overlap layers for AS preconditioner
HALO ! AS restriction operator: NONE HALO
NONE ! AS prolongation operator: NONE SUM AVG
@ -32,7 +32,7 @@ LLK ! AINV variant
1 ! Inverse Fill level P for INVK
1.d-4 ! Threshold T for ILU(T,P)
%%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%%
L1-JACOBI ! Second (post) smoother, ignored if NONE
NONE ! Second (post) smoother, ignored if NONE
6 ! Number of sweeps for (post) smoother
1 ! degree for polynomial smoother
POLY_LOTTES_BETA ! Polynomial variant

Loading…
Cancel
Save