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

Adjustments for POLYNOMIAL smoothers.

Browse Source
This commit is contained in:
sfilippone
2023-11-17 13:14:06 +01:00
parent 2dd1cbd3dc
commit ec9fcb1bcc
22 changed files with 200 additions and 81 deletions
Download Patch File Download Diff File Expand all files Collapse all files
amgprec/Makefile
+2 -2
View File
@@ -9,7 +9,7 @@ FINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(INCDIR) $(PSBLAS_INCLUDES)
DMODOBJS=amg_d_prec_type.o \ DMODOBJS=amg_d_prec_type.o \
amg_d_inner_mod.o amg_d_ilu_solver.o amg_d_diag_solver.o amg_d_jac_smoother.o amg_d_as_smoother.o \ amg_d_inner_mod.o amg_d_ilu_solver.o amg_d_diag_solver.o amg_d_jac_smoother.o amg_d_as_smoother.o \
amg_d_poly_smoother.o amg_d_beta_coeff_mod.o\ amg_d_poly_smoother.o amg_d_poly_coeff_mod.o\
amg_d_umf_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o amg_d_id_solver.o\ amg_d_umf_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o amg_d_id_solver.o\
amg_d_base_solver_mod.o amg_d_base_smoother_mod.o amg_d_onelev_mod.o \ amg_d_base_solver_mod.o amg_d_base_smoother_mod.o amg_d_onelev_mod.o \
amg_d_gs_solver.o amg_d_mumps_solver.o amg_d_jac_solver.o \ amg_d_gs_solver.o amg_d_mumps_solver.o amg_d_jac_solver.o \
@@ -165,7 +165,7 @@ amg_d_jac_smoother.o: amg_d_diag_solver.o
amg_dprecinit.o amg_dprecset.o: amg_d_diag_solver.o amg_d_ilu_solver.o \ 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_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_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_beta_coeff_mod.o amg_d_poly_smoother.o: amg_d_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_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 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
amgprec/amg_base_prec_type.F90
+12
View File
@@ -321,6 +321,12 @@ module amg_base_prec_type
integer(psb_ipk_), parameter :: amg_repl_mat_ = 1 integer(psb_ipk_), parameter :: amg_repl_mat_ = 1
integer(psb_ipk_), parameter :: amg_max_coarse_mat_ = amg_repl_mat_ integer(psb_ipk_), parameter :: amg_max_coarse_mat_ = amg_repl_mat_
! !
! Legal values for entry: amg_poly_variant_
!
integer(psb_ipk_), parameter :: amg_poly_lottes_ = 0
integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1
integer(psb_ipk_), parameter :: amg_poly_new_ = 2
!
! Legal values for entry: amg_prec_status_ ! Legal values for entry: amg_prec_status_
! !
integer(psb_ipk_), parameter :: amg_prec_built_ = 98765 integer(psb_ipk_), parameter :: amg_prec_built_ = 98765
@@ -560,6 +566,12 @@ contains
val = amg_as_ val = amg_as_
case('POLY') case('POLY')
val = amg_poly_ val = amg_poly_
case('POLY_LOTTES')
val = amg_poly_lottes_
case('POLY_LOTTES_BETA')
val = amg_poly_lottes_beta_
case('POLY_NEW')
val = amg_poly_new_
case('A_NORMI') case('A_NORMI')
val = amg_max_norm_ val = amg_max_norm_
case('USER_CHOICE') case('USER_CHOICE')
amgprec/amg_d_poly_coeff_mod.f90
+2 -2
View File
@@ -49,7 +49,7 @@
! main diagonal block), so that it becomes possible to implement ! main diagonal block), so that it becomes possible to implement
! a pure Jacobi or L1-Jacobi global solver. ! a pure Jacobi or L1-Jacobi global solver.
! !
module amg_d_beta_coeff_mod module amg_d_poly_coeff_mod
use psb_base_mod use psb_base_mod
real(psb_dpk_), parameter :: amg_d_beta_vect(900) = [ & real(psb_dpk_), parameter :: amg_d_beta_vect(900) = [ &
@@ -513,4 +513,4 @@ module amg_d_beta_coeff_mod
real(psb_dpk_), parameter :: amg_d_beta_mat(30,30)=reshape(amg_d_beta_vect,[30,30]) real(psb_dpk_), parameter :: amg_d_beta_mat(30,30)=reshape(amg_d_beta_vect,[30,30])
end module amg_d_beta_coeff_mod end module amg_d_poly_coeff_mod
amgprec/amg_d_poly_smoother.f90
+3 -2
View File
@@ -51,14 +51,14 @@
! !
module amg_d_poly_smoother module amg_d_poly_smoother
use amg_d_base_smoother_mod use amg_d_base_smoother_mod
use amg_d_beta_coeff_mod use amg_d_poly_coeff_mod
type, extends(amg_d_base_smoother_type) :: amg_d_poly_smoother_type type, extends(amg_d_base_smoother_type) :: amg_d_poly_smoother_type
! The local solver component is inherited from the ! The local solver component is inherited from the
! parent type. ! parent type.
! class(amg_d_base_solver_type), allocatable :: sv ! class(amg_d_base_solver_type), allocatable :: sv
! !
integer(psb_ipk_) :: pdegree integer(psb_ipk_) :: pdegree, variant
type(psb_dspmat_type), pointer :: pa => null() type(psb_dspmat_type), pointer :: pa => null()
real(psb_dpk_), allocatable :: poly_beta(:) real(psb_dpk_), allocatable :: poly_beta(:)
real(psb_dpk_) :: rho_ba real(psb_dpk_) :: rho_ba
@@ -319,6 +319,7 @@ contains
! !
sm%pdegree = 1 sm%pdegree = 1
sm%rho_ba = dzero sm%rho_ba = dzero
sm%variant = amg_poly_lottes_
if (allocated(sm%sv)) then if (allocated(sm%sv)) then
call sm%sv%default() call sm%sv%default()
amgprec/amg_d_prec_mod.f90
-1
View File
@@ -47,7 +47,6 @@ module amg_d_prec_mod
use amg_d_prec_type use amg_d_prec_type
use amg_d_jac_smoother use amg_d_jac_smoother
use amg_d_as_smoother use amg_d_as_smoother
use amg_d_poly_smoother
use amg_d_id_solver use amg_d_id_solver
use amg_d_diag_solver use amg_d_diag_solver
use amg_d_l1_diag_solver use amg_d_l1_diag_solver
amgprec/impl/amg_dcprecset.F90
+1 -7
View File
@@ -81,7 +81,6 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx)
use amg_d_prec_mod, amg_protect_name => amg_dcprecseti use amg_d_prec_mod, amg_protect_name => amg_dcprecseti
use amg_d_jac_smoother use amg_d_jac_smoother
use amg_d_as_smoother use amg_d_as_smoother
use amg_d_poly_smoother
use amg_d_diag_solver use amg_d_diag_solver
use amg_d_l1_diag_solver use amg_d_l1_diag_solver
use amg_d_ilu_solver use amg_d_ilu_solver
@@ -126,7 +125,7 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx)
info = 3111 info = 3111
write(psb_err_unit,*) name,& write(psb_err_unit,*) name,&
& ': Error: uninitialized preconditioner,',& & ': Error: uninitialized preconditioner,',&
& ' should call amg_PRECINIT' &' should call amg_PRECINIT'
return return
endif endif
@@ -313,7 +312,6 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
use amg_d_prec_mod, amg_protect_name => amg_dcprecsetc use amg_d_prec_mod, amg_protect_name => amg_dcprecsetc
use amg_d_jac_smoother use amg_d_jac_smoother
use amg_d_as_smoother use amg_d_as_smoother
use amg_d_poly_smoother
use amg_d_diag_solver use amg_d_diag_solver
use amg_d_l1_diag_solver use amg_d_l1_diag_solver
use amg_d_ilu_solver use amg_d_ilu_solver
@@ -404,10 +402,6 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
do il=ilev_, ilmax_ do il=ilev_, ilmax_
call p%precv(il)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) call p%precv(il)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
end do end do
case ('POLY')
do il=ilev_, ilmax_
call p%precv(il)%set('SMOOTHER_TYPE',amg_poly_,info,pos=pos)
end do
case ('L1-BJAC') case ('L1-BJAC')
do il=ilev_, ilmax_ do il=ilev_, ilmax_
call p%precv(il)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos) call p%precv(il)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
amgprec/impl/amg_dprecinit.F90
+1 -1
View File
@@ -163,7 +163,7 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
allocate(prec%precv(nlev_),stat=info) allocate(prec%precv(nlev_),stat=info)
allocate(amg_d_poly_smoother_type :: prec%precv(ilev_)%sm, stat=info) allocate(amg_d_poly_smoother_type :: prec%precv(ilev_)%sm, stat=info)
if (info /= psb_success_) return if (info /= psb_success_) return
allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) allocate(amg_d_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
call prec%precv(ilev_)%default() call prec%precv(ilev_)%default()
case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
amgprec/impl/level/amg_d_base_onelev_csetc.F90
+2 -7
View File
@@ -207,16 +207,11 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
case ('NONE','NOPREC','FACT_NONE') case ('NONE','NOPREC','FACT_NONE')
call lv%set(amg_d_id_solver_mold,info,pos=pos) call lv%set(amg_d_id_solver_mold,info,pos=pos)
case ('DIAG') case ('DIAG','JACOBI')
call lv%set(amg_d_diag_solver_mold,info,pos=pos) call lv%set(amg_d_diag_solver_mold,info,pos=pos)
case ('JACOBI') case ('L1-DIAG','L1-JACOBI')
call lv%set(amg_d_jac_solver_mold,info,pos=pos)
case ('L1-DIAG')
call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos) call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos)
case ('L1-JACOBI')
call lv%set(amg_d_l1_jac_solver_mold,info,pos=pos)
case ('GS','FGS','FWGS') case ('GS','FGS','FWGS')
call lv%set(amg_d_gs_solver_mold,info,pos=pos) call lv%set(amg_d_gs_solver_mold,info,pos=pos)
amgprec/impl/level/amg_d_base_onelev_cseti.F90
-6
View File
@@ -45,7 +45,6 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx)
use amg_d_parmatch_aggregator_mod use amg_d_parmatch_aggregator_mod
use amg_d_jac_smoother use amg_d_jac_smoother
use amg_d_as_smoother use amg_d_as_smoother
use amg_d_poly_smoother
use amg_d_diag_solver use amg_d_diag_solver
use amg_d_l1_diag_solver use amg_d_l1_diag_solver
use amg_d_ilu_solver use amg_d_ilu_solver
@@ -80,7 +79,6 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx)
type(amg_d_jac_smoother_type) :: amg_d_jac_smoother_mold 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_l1_jac_smoother_type) :: amg_d_l1_jac_smoother_mold
type(amg_d_as_smoother_type) :: amg_d_as_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_diag_solver_type) :: amg_d_diag_solver_mold
type(amg_d_l1_diag_solver_type) :: amg_d_l1_diag_solver_mold type(amg_d_l1_diag_solver_type) :: amg_d_l1_diag_solver_mold
type(amg_d_ilu_solver_type) :: amg_d_ilu_solver_mold type(amg_d_ilu_solver_type) :: amg_d_ilu_solver_mold
@@ -143,10 +141,6 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx)
call lv%set(amg_d_as_smoother_mold,info,pos=pos) call lv%set(amg_d_as_smoother_mold,info,pos=pos)
if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos) if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos)
case (amg_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 (amg_fbgs_) case (amg_fbgs_)
call lv%set(amg_d_jac_smoother_mold,info,pos='pre') 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') if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre')
amgprec/impl/smoother/amg_d_jac_smoother_apply_vect.f90
+10 -3
View File
@@ -109,7 +109,7 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
endif endif
if(sm%checkres) then if(.true..or.sm%checkres) then
call psb_geall(r,desc_data,info) call psb_geall(r,desc_data,info)
call psb_geasb(r,desc_data,info) call psb_geasb(r,desc_data,info)
resdenum = psb_genrm2(x,desc_data,info) resdenum = psb_genrm2(x,desc_data,info)
@@ -159,7 +159,10 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& a_err='wrong init to smoother_apply') & a_err='wrong init to smoother_apply')
goto 9999 goto 9999
end select 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 do i=1, sweeps-1
! !
! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)),
@@ -173,9 +176,13 @@ 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') call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y')
if (info /= psb_success_) exit 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 if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then
call psb_geaxpby(done,x,dzero,r,r,desc_data,info) call psb_geaxpby(done,x,dzero,r,desc_data,info)
call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info)
res = psb_genrm2(r,desc_data,info) res = psb_genrm2(r,desc_data,info)
if( sm%printres ) then if( sm%printres ) then
amgprec/impl/smoother/amg_d_jac_smoother_clone.f90
-1
View File
@@ -72,7 +72,6 @@ subroutine amg_d_jac_smoother_clone(sm,smout,info)
smo%checkiter = sm%checkiter smo%checkiter = sm%checkiter
smo%printiter = sm%printiter smo%printiter = sm%printiter
smo%tol = sm%tol smo%tol = sm%tol
smo%pa => sm%pa
call sm%nd%clone(smo%nd,info) call sm%nd%clone(smo%nd,info)
if ((info==psb_success_).and.(allocated(sm%sv))) then if ((info==psb_success_).and.(allocated(sm%sv))) then
allocate(smout%sv,mold=sm%sv,stat=info) allocate(smout%sv,mold=sm%sv,stat=info)
amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90
+64 -15
View File
@@ -129,15 +129,64 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
call ty%zero() call ty%zero()
call tz%zero() call tz%zero()
do i=1, sm%pdegree select case(sm%variant)
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Y') case(amg_poly_lottes_)
cz = (2*i*done-3)/(2*i*done+done) ! b == x
cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) ! x == tx
call psb_geaxpby(cr,ty,cz,tz,desc_data,info) !
call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info) do i=1, sm%pdegree
call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_) ! B r_{k-1}
end do 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)
cr = (8*i*done-4)/((2*i*done+done)*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(done,tz,done,tx,desc_data,info)
if (.false.) then
call psb_geaxpby(done,x,dzero,r,desc_data,info)
call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
else
call psb_spmm(-done,sm%pa,tz,done,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
case(amg_poly_lottes_beta_)
! b == x
! x == tx
!
do i=1, sm%pdegree
! 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)
cr = (8*i*done-4)/((2*i*done+done)*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,done,tx,desc_data,info)
if (.false.) then
call psb_geaxpby(done,x,dzero,r,desc_data,info)
call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
else
call psb_spmm(-done,sm%pa,tz,done,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
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_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@@ -339,19 +388,19 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
!!$ !!$
!!$ endif !!$ endif
!!$ !!$
if (.not.(4*n_col <= size(work))) then if (.not.(4*n_col <= size(work))) then
deallocate(aux) deallocate(aux)
endif endif
!!$ if(sm%checkres) then !!$ if(sm%checkres) then
!!$ call psb_gefree(r,desc_data,info) !!$ call psb_gefree(r,desc_data,info)
!!$ end if !!$ end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine amg_d_poly_smoother_apply_vect end subroutine amg_d_poly_smoother_apply_vect
amgprec/impl/smoother/amg_d_poly_smoother_bld.f90
+45 -4
View File
@@ -39,7 +39,8 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
use psb_base_mod use psb_base_mod
use amg_d_diag_solver use amg_d_diag_solver
use amg_d_beta_coeff_mod use amg_d_l1_diag_solver
use amg_d_poly_coeff_mod
use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_bld use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_bld
Implicit None Implicit None
@@ -55,6 +56,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
type(psb_dspmat_type) :: tmpa type(psb_dspmat_type) :: tmpa
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros
type(psb_ctxt_type) :: ctxt type(psb_ctxt_type) :: ctxt
real(psb_dpk_), allocatable :: da(:), dsv(:)
integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level
character(len=20) :: name='d_poly_smoother_bld', ch_err character(len=20) :: name='d_poly_smoother_bld', ch_err
@@ -83,15 +85,54 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
goto 9999 goto 9999
end if end if
sm%pa => a 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) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='sv%build') & a_err='sv%build')
goto 9999 goto 9999
end if end if
if (sm%rho_ba <= dzero) then if (.true.) then
sm%rho_ba = psb_dspnrm1(a,desc_a,info) select type(ssv => sm%sv)
class is(amg_d_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 = done
end select
else
block
type(psb_d_vect_type) :: tq, tz,wv(2)
real(psb_dpk_) :: qnrm, 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(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
end do
sm%rho_ba = lambda
end block
end if end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90
+4 -2
View File
@@ -53,8 +53,10 @@ subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select case(psb_toupper(trim(what))) select case(psb_toupper(trim(what)))
case default case('POLY-VARIANT')
call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx) 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)
end select end select
if (info /= psb_success_) then if (info /= psb_success_) then
amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90
+9 -1
View File
@@ -54,8 +54,16 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
select case(psb_toupper(what)) select case(psb_toupper(what))
case('SMOOTHER_DEGREE') case('POLY_DEGREE')
sm%pdegree = val 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_'
sm%variant = amg_poly_lottes_
end select
case default case default
call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx) call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx)
end select end select
amgprec/impl/smoother/amg_d_poly_smoother_descr.f90
+8
View File
@@ -78,6 +78,14 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix)
end if end if
write(iout_,*) trim(prefix_), ' Polynomial smoother ' write(iout_,*) trim(prefix_), ' Polynomial smoother '
select case(sm%variant)
case(amg_poly_lottes_)
write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES'
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_), ' Degree: ',sm%pdegree
write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba 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)
amgprec/impl/solver/amg_c_ilu_solver_bld.f90
+3 -3
View File
@@ -99,11 +99,11 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
select case(sv%fact_type) select case(sv%fact_type)
case (amg_ilu_n_) case (amg_ilu_n_)
psb_fctype = amg_ilu_n_ psb_fctype = psb_ilu_n_
case (amg_milu_n_) case (amg_milu_n_)
psb_fctype = amg_milu_n_ psb_fctype = psb_milu_n_
case (amg_ilu_t_) case (amg_ilu_t_)
psb_fctype = amg_ilu_t_ psb_fctype = psb_ilu_t_
case default case default
! If we end up here, something was wrong up in the call chain. ! If we end up here, something was wrong up in the call chain.
info = psb_err_input_value_invalid_i_ info = psb_err_input_value_invalid_i_
amgprec/impl/solver/amg_d_ilu_solver_bld.f90
+3 -3
View File
@@ -99,11 +99,11 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
select case(sv%fact_type) select case(sv%fact_type)
case (amg_ilu_n_) case (amg_ilu_n_)
psb_fctype = amg_ilu_n_ psb_fctype = psb_ilu_n_
case (amg_milu_n_) case (amg_milu_n_)
psb_fctype = amg_milu_n_ psb_fctype = psb_milu_n_
case (amg_ilu_t_) case (amg_ilu_t_)
psb_fctype = amg_ilu_t_ psb_fctype = psb_ilu_t_
case default case default
! If we end up here, something was wrong up in the call chain. ! If we end up here, something was wrong up in the call chain.
info = psb_err_input_value_invalid_i_ info = psb_err_input_value_invalid_i_
amgprec/impl/solver/amg_s_ilu_solver_bld.f90
+3 -3
View File
@@ -99,11 +99,11 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
select case(sv%fact_type) select case(sv%fact_type)
case (amg_ilu_n_) case (amg_ilu_n_)
psb_fctype = amg_ilu_n_ psb_fctype = psb_ilu_n_
case (amg_milu_n_) case (amg_milu_n_)
psb_fctype = amg_milu_n_ psb_fctype = psb_milu_n_
case (amg_ilu_t_) case (amg_ilu_t_)
psb_fctype = amg_ilu_t_ psb_fctype = psb_ilu_t_
case default case default
! If we end up here, something was wrong up in the call chain. ! If we end up here, something was wrong up in the call chain.
info = psb_err_input_value_invalid_i_ info = psb_err_input_value_invalid_i_
amgprec/impl/solver/amg_z_ilu_solver_bld.f90
+3 -3
View File
@@ -99,11 +99,11 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
select case(sv%fact_type) select case(sv%fact_type)
case (amg_ilu_n_) case (amg_ilu_n_)
psb_fctype = amg_ilu_n_ psb_fctype = psb_ilu_n_
case (amg_milu_n_) case (amg_milu_n_)
psb_fctype = amg_milu_n_ psb_fctype = psb_milu_n_
case (amg_ilu_t_) case (amg_ilu_t_)
psb_fctype = amg_ilu_t_ psb_fctype = psb_ilu_t_
case default case default
! If we end up here, something was wrong up in the call chain. ! If we end up here, something was wrong up in the call chain.
info = psb_err_input_value_invalid_i_ info = psb_err_input_value_invalid_i_
samples/advanced/pdegen/amg_d_pde3d.F90
+17 -9
View File
@@ -145,6 +145,7 @@ program amg_d_pde3d
character(len=16) :: smther ! (pre-)smoother type: BJAC, AS character(len=16) :: smther ! (pre-)smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps
integer(psb_ipk_) :: degree ! degree for polynomial smoother integer(psb_ipk_) :: degree ! degree for polynomial smoother
character(len=16) :: pvariant ! Polynomial smoother variant
integer(psb_ipk_) :: novr ! number of overlap layers integer(psb_ipk_) :: novr ! number of overlap layers
character(len=16) :: restr ! restriction over application of AS character(len=16) :: restr ! restriction over application of AS
character(len=16) :: prol ! prolongation over application of AS character(len=16) :: prol ! prolongation over application of AS
@@ -154,12 +155,13 @@ program amg_d_pde3d
character(len=16) :: variant ! AINV variant: LLK, etc character(len=16) :: variant ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill ! Inverse fill-in for INVK integer(psb_ipk_) :: invfill ! Inverse fill-in for INVK
real(psb_dpk_) :: thr ! threshold for ILUT factorization real(psb_dpk_) :: thr ! threshold for ILUT factorization
! AMG post-smoother; ignored by 1-lev preconditioner ! AMG post-smoother; ignored by 1-lev preconditioner
character(len=16) :: smther2 ! post-smoother type: BJAC, AS character(len=16) :: smther2 ! post-smoother type: BJAC, AS
integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps
integer(psb_ipk_) :: degree2 ! degree for polynomial smoother integer(psb_ipk_) :: degree2 ! degree for polynomial smoother
character(len=16) :: pvariant2 ! Polynomial smoother variant
integer(psb_ipk_) :: novr2 ! number of overlap layers integer(psb_ipk_) :: novr2 ! number of overlap layers
character(len=16) :: restr2 ! restriction over application of AS character(len=16) :: restr2 ! restriction over application of AS
character(len=16) :: prol2 ! prolongation over application of AS character(len=16) :: prol2 ! prolongation over application of AS
@@ -169,7 +171,7 @@ program amg_d_pde3d
character(len=16) :: variant2 ! AINV variant: LLK, etc character(len=16) :: variant2 ! AINV variant: LLK, etc
integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization
integer(psb_ipk_) :: invfill2 ! Inverse fill-in for INVK integer(psb_ipk_) :: invfill2 ! Inverse fill-in for INVK
real(psb_dpk_) :: thr2 ! threshold for ILUT factorization real(psb_dpk_) :: thr2 ! threshold for ILUT factorization
! coarsest-level solver ! coarsest-level solver
character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST
@@ -291,7 +293,8 @@ program amg_d_pde3d
call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('sub_solve', p_choice%solve, info) call prec%set('sub_solve', p_choice%solve, info)
call prec%set('solver_sweeps', p_choice%ssweeps, info) call prec%set('solver_sweeps', p_choice%ssweeps, info)
call prec%set('smoother_degree', p_choice%degree, info) call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%variant, info)
if (psb_toupper(p_choice%solve)=='MUMPS') & if (psb_toupper(p_choice%solve)=='MUMPS') &
& call prec%set('mumps_loc_glob','local_solver',info) & call prec%set('mumps_loc_glob','local_solver',info)
call prec%set('sub_fillin', p_choice%fill, info) call prec%set('sub_fillin', p_choice%fill, info)
@@ -339,7 +342,8 @@ program amg_d_pde3d
call prec%set('smoother_type', p_choice%smther, info) call prec%set('smoother_type', p_choice%smther, info)
call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info)
call prec%set('smoother_degree', p_choice%degree, info) call prec%set('poly_degree', p_choice%degree, info)
call prec%set('poly_variant', p_choice%variant, info)
select case (psb_toupper(p_choice%smther)) select case (psb_toupper(p_choice%smther))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
@@ -370,7 +374,8 @@ program amg_d_pde3d
if (psb_toupper(p_choice%smther2) /= 'NONE') then if (psb_toupper(p_choice%smther2) /= 'NONE') then
call prec%set('smoother_type', p_choice%smther2, info,pos='post') 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('smoother_sweeps', p_choice%jsweeps2, info,pos='post')
call prec%set('smoother_degree', p_choice%degree2, 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')
select case (psb_toupper(p_choice%smther2)) select case (psb_toupper(p_choice%smther2))
case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
! do nothing ! do nothing
@@ -407,7 +412,6 @@ program amg_d_pde3d
call prec%set('coarse_sweeps', p_choice%cjswp, info) call prec%set('coarse_sweeps', p_choice%cjswp, info)
end select end select
!!$ call prec%descr(info,iout=psb_out_unit)
! build the preconditioner ! build the preconditioner
call psb_barrier(ctxt) call psb_barrier(ctxt)
@@ -587,7 +591,8 @@ contains
! First smoother / 1-lev preconditioner ! First smoother / 1-lev preconditioner
call read_data(prec%smther,inp_unit) ! smoother type 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%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%degree,inp_unit) ! Degree of Polynomial smoother
call read_data(prec%variant,inp_unit) ! variant for Polynomial
call read_data(prec%novr,inp_unit) ! number of overlap layers 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%restr,inp_unit) ! restriction over application of AS
call read_data(prec%prol,inp_unit) ! prolongation over application of AS call read_data(prec%prol,inp_unit) ! prolongation over application of AS
@@ -600,12 +605,13 @@ contains
! Second smoother/ AMG post-smoother (if NONE ignored in main) ! Second smoother/ AMG post-smoother (if NONE ignored in main)
call read_data(prec%smther2,inp_unit) ! smoother type call read_data(prec%smther2,inp_unit) ! smoother type
call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps
call read_data(prec%degree2,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%novr2,inp_unit) ! number of overlap layers 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%restr2,inp_unit) ! restriction over application of AS
call read_data(prec%prol2,inp_unit) ! prolongation over application of AS call read_data(prec%prol2,inp_unit) ! prolongation over application of AS
call read_data(prec%solve2,inp_unit) ! local subsolver call read_data(prec%solve2,inp_unit) ! local subsolver
call read_data(prec%ssweeps2,inp_unit) ! inner solver sweeps call read_data(prec%ssweeps2,inp_unit) ! inner solver sweeps
call read_data(prec%variant2,inp_unit) ! AINV variant call read_data(prec%variant2,inp_unit) ! AINV variant
call read_data(prec%fill2,inp_unit) ! fill-in for incomplete LU call read_data(prec%fill2,inp_unit) ! fill-in for incomplete LU
call read_data(prec%invfill2,inp_unit) !Inverse fill-in for INVK call read_data(prec%invfill2,inp_unit) !Inverse fill-in for INVK
@@ -672,6 +678,7 @@ contains
call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%smther)
call psb_bcast(ctxt,prec%jsweeps) call psb_bcast(ctxt,prec%jsweeps)
call psb_bcast(ctxt,prec%degree) call psb_bcast(ctxt,prec%degree)
call psb_bcast(ctxt,prec%variant)
call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%novr)
call psb_bcast(ctxt,prec%restr) call psb_bcast(ctxt,prec%restr)
call psb_bcast(ctxt,prec%prol) call psb_bcast(ctxt,prec%prol)
@@ -685,6 +692,7 @@ contains
call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%smther2)
call psb_bcast(ctxt,prec%jsweeps2) call psb_bcast(ctxt,prec%jsweeps2)
call psb_bcast(ctxt,prec%degree2) call psb_bcast(ctxt,prec%degree2)
call psb_bcast(ctxt,prec%variant2)
call psb_bcast(ctxt,prec%novr2) call psb_bcast(ctxt,prec%novr2)
call psb_bcast(ctxt,prec%restr2) call psb_bcast(ctxt,prec%restr2)
call psb_bcast(ctxt,prec%prol2) call psb_bcast(ctxt,prec%prol2)
samples/advanced/pdegen/runs/amg_pde3d.inp
+8 -6
View File
@@ -4,16 +4,17 @@ CSR ! Storage format CSR COO JAD
CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE
BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
2 ! ISTOPC 2 ! ISTOPC
00500 ! ITMAX 00050 ! ITMAX
1 ! ITRACE 1 ! ITRACE
30 ! IRST (restart for RGMRES and BiCGSTABL) 30 ! IRST (restart for RGMRES and BiCGSTABL)
1.d-6 ! EPS 1.d-6 ! EPS
%%%%%%%%%%% Main preconditioner choices %%%%%%%%%%%%%%%% %%%%%%%%%%% Main preconditioner choices %%%%%%%%%%%%%%%%
ML-VBM-VCYCLE-FBGS-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars) ML-VBM-VCYCLE-FBGS-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars)
BJAC ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML
%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% %%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
BJAC ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous.
1 ! Number of sweeps for smoother 4 ! Number of sweeps for smoother
4 ! degree for polynomial smoother
0 ! Number of overlap layers for AS preconditioner 0 ! Number of overlap layers for AS preconditioner
HALO ! AS restriction operator: NONE HALO HALO ! AS restriction operator: NONE HALO
NONE ! AS prolongation operator: NONE SUM AVG NONE ! AS prolongation operator: NONE SUM AVG
@@ -24,8 +25,9 @@ LLK ! AINV variant
1 ! Inverse Fill level P for INVK 1 ! Inverse Fill level P for INVK
1.d-4 ! Threshold T for ILU(T,P) 1.d-4 ! Threshold T for ILU(T,P)
%%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%% %%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%%
NONE ! Second (post) smoother, ignored if NONE NONE ! Second (post) smoother, ignored if NONE
1 ! Number of sweeps for (post) smoother 4 ! Number of sweeps for (post) smoother
4 ! degree for polynomial smoother
0 ! Number of overlap layers for AS preconditioner 0 ! Number of overlap layers for AS preconditioner
HALO ! AS restriction operator: NONE HALO HALO ! AS restriction operator: NONE HALO
NONE ! AS prolongation operator: NONE SUM AVG NONE ! AS prolongation operator: NONE SUM AVG