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

Temporary checkpoint, working version, to be investigated further.

Browse Source
This commit is contained in:
sfilippone
2023-11-17 14:45:18 +01:00
parent 14cd4cde76
commit bb262275a1
7 changed files with 63 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files
amgprec/amg_base_prec_type.F90
+1 -1
View File
@@ -464,7 +464,7 @@ contains
character(len=*), parameter :: name='amg_stringval'
! Local variable
integer :: index_tab
character(len=15) ::string2
character(len=128) ::string2
index_tab=index(string,char(9))
if (index_tab.NE.0) then
string2=string(1:index_tab-1)
amgprec/impl/smoother/amg_d_poly_smoother_bld.f90
+14 -11
View File
@@ -98,7 +98,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
goto 9999
end if
if (.true.) then
if (.false.) then
select type(ssv => sm%sv)
class is(amg_d_l1_diag_solver_type)
da = a%arwsum(info)
@@ -110,28 +110,31 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
end select
else
block
type(psb_d_vect_type) :: tq, tz,wv(2)
real(psb_dpk_) :: qnrm, lambda
type(psb_d_vect_type) :: tq, tt, tz,wv(2)
real(psb_dpk_) :: znrm, lambda
real(psb_dpk_),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(done)
call psb_geasb(tq,desc_a,info,mold=vmold)
call psb_spmm(done,a,tq,dzero,tz,desc_a,info) ! z_1 = A q_0
do i=1,10
call sm%sv%apply_v(done,tz,dzero,tq,desc_a,'NoTrans',work,wv,info) ! q_k = M^{-1} q_k
qnrm = psb_genrmi(tq,desc_a,info) ! qnrm = |q_k|_inf
call tq%scal((done/qnrm)) ! q_k = q_k/qnrm
call psb_spmm(done,a,tq,dzero,tz,desc_a,info) ! z_{k=1} = A q_k
lambda = psb_gedot(tq,tz,desc_a,info) ! lambda = q_k^T z_{k+1} = q_k^T A q_k
write(0,*) 'BLD: lambda estimate ',i,lambda
call psb_spmm(done,a,tq,dzero,tt,desc_a,info) !
call sm%sv%apply_v(done,tt,dzero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = BA q_k
do i=1,20
znrm = psb_genrm2(tz,desc_a,info) ! znrm = |z_k|_2
call psb_geaxpby((done/znrm),tz,dzero,tq,desc_a,info) ! q_k = z_k/znrm
call psb_spmm(done,a,tq,dzero,tt,desc_a,info) ! t_{k+1} = BA q_k
call sm%sv%apply_v(done,tt,dzero,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
sm%rho_ba = done
end block
end if
amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90
+1 -1
View File
@@ -53,7 +53,7 @@ subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx)
call psb_erractionsave(err_act)
select case(psb_toupper(trim(what)))
case('POLY-VARIANT')
case('POLY_VARIANT')
call sm%set(what,amg_stringval(val),info,idx=idx)
case default
call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx)
amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90
+1 -1
View File
@@ -61,7 +61,7 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx)
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_'
write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val
sm%variant = amg_poly_lottes_
end select
case default
amgprec/impl/smoother/amg_d_poly_smoother_descr.f90
+5 -3
View File
@@ -81,14 +81,16 @@ 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_), ' rho_ba: ',sm%rho_ba
case(amg_poly_lottes_beta_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES_BETA'
case default
write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???'
end select
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 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)
samples/advanced/pdegen/amg_d_pde3d.F90
+32 -31
View File
@@ -80,7 +80,7 @@ program amg_d_pde3d
implicit none
! input parameters
character(len=20) :: kmethd, ptype
character(len=24) :: kmethd, ptype
character(len=5) :: afmt, pdecoeff
integer(psb_ipk_) :: idim
integer(psb_epk_) :: system_size
@@ -120,21 +120,21 @@ program amg_d_pde3d
! preconditioner type
character(len=40) :: descr ! verbose description of the prec
character(len=10) :: ptype ! preconditioner type
character(len=24) :: ptype ! preconditioner type
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=24) :: 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=24) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED
character(len=24) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC
character(len=24) :: 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=24) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE
character(len=24) :: 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
@@ -142,43 +142,43 @@ program amg_d_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=24) :: 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=16) :: pvariant ! Polynomial smoother variant
character(len=24) :: pvariant ! Polynomial smoother 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=24) :: restr ! restriction over application of AS
character(len=24) :: prol ! prolongation over application of AS
character(len=24) :: 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=24) :: 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=24) :: smther2 ! post-smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps
integer(psb_ipk_) :: degree2 ! degree for polynomial smoother
character(len=16) :: pvariant2 ! Polynomial smoother variant
character(len=24) :: pvariant2 ! Polynomial smoother 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=24) :: restr2 ! restriction over application of AS
character(len=24) :: prol2 ! prolongation over application of AS
character(len=24) :: 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=24) :: 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=24) :: cmat ! coarsest matrix layout: REPL, DIST
character(len=24) :: 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=24) :: 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
@@ -200,7 +200,7 @@ program amg_d_pde3d
! other variables
integer(psb_ipk_) :: info, i, k
character(len=20) :: name,ch_err
character(len=24) :: name,ch_err
type(psb_d_csr_sparse_mat) :: amold
info=psb_success_
@@ -294,7 +294,7 @@ program amg_d_pde3d
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%variant, 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)
@@ -343,7 +343,8 @@ program amg_d_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%variant, info)
write(0,*) 'pvariant :',p_choice%pvariant
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')
@@ -375,7 +376,7 @@ program amg_d_pde3d
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%variant2, 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
@@ -592,7 +593,7 @@ contains
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) ! Degree of Polynomial smoother
call read_data(prec%variant,inp_unit) ! variant for Polynomial
call read_data(prec%pvariant,inp_unit) ! variant for Polynomial
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
@@ -606,7 +607,7 @@ contains
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) ! Degree of Polynomial smoother
call read_data(prec%variant2,inp_unit) ! Polynomial smoother variant
call read_data(prec%pvariant2,inp_unit) ! Polynomial smoother variant
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
@@ -678,7 +679,7 @@ contains
call psb_bcast(ctxt,prec%smther)
call psb_bcast(ctxt,prec%jsweeps)
call psb_bcast(ctxt,prec%degree)
call psb_bcast(ctxt,prec%variant)
call psb_bcast(ctxt,prec%pvariant)
call psb_bcast(ctxt,prec%novr)
call psb_bcast(ctxt,prec%restr)
call psb_bcast(ctxt,prec%prol)
@@ -692,7 +693,7 @@ contains
call psb_bcast(ctxt,prec%smther2)
call psb_bcast(ctxt,prec%jsweeps2)
call psb_bcast(ctxt,prec%degree2)
call psb_bcast(ctxt,prec%variant2)
call psb_bcast(ctxt,prec%pvariant2)
call psb_bcast(ctxt,prec%novr2)
call psb_bcast(ctxt,prec%restr2)
call psb_bcast(ctxt,prec%prol2)
samples/advanced/pdegen/runs/amg_pde3d.inp
+5 -5
View File
@@ -1,6 +1,6 @@
%%%%%%%%%%% General arguments % Lines starting with % are ignored.
CSR ! Storage format CSR COO JAD
0050 ! IDIM; domain size. Linear system size is IDIM**3
0150 ! IDIM; domain size. Linear system size is IDIM**3
CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE
BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
2 ! ISTOPC
@@ -12,8 +12,8 @@ BICGSTAB ! 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
%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous.
1 ! Number of sweeps for smoother
L1-JACOBI ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 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
@@ -26,8 +26,8 @@ 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 %%%%%%%%%%%%%%%%
FBGS ! Second (post) smoother, ignored if NONE
1 ! Number of sweeps for (post) smoother
L1-JACOBI ! Second (post) smoother, ignored if NONE
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