From 79317cb392e473905d7bb363d28dfc63db848f64 Mon Sep 17 00:00:00 2001 From: sfilippone Date: Fri, 17 Nov 2023 15:35:15 +0100 Subject: [PATCH] Additional fields for rho(BA) estimate. --- amgprec/amg_base_prec_type.F90 | 5 ++ amgprec/amg_d_poly_smoother.f90 | 13 +-- .../impl/smoother/amg_d_poly_smoother_bld.f90 | 80 ++++++++++--------- .../smoother/amg_d_poly_smoother_csetc.f90 | 2 + .../smoother/amg_d_poly_smoother_cseti.f90 | 15 ++++ .../smoother/amg_d_poly_smoother_csetr.f90 | 9 ++- samples/advanced/pdegen/runs/amg_pde3d.inp | 8 +- 7 files changed, 87 insertions(+), 45 deletions(-) diff --git a/amgprec/amg_base_prec_type.F90 b/amgprec/amg_base_prec_type.F90 index 59b1acce..5d75e274 100644 --- a/amgprec/amg_base_prec_type.F90 +++ b/amgprec/amg_base_prec_type.F90 @@ -326,6 +326,9 @@ module amg_base_prec_type 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 + + integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0 + ! ! Legal values for entry: amg_prec_status_ ! @@ -572,6 +575,8 @@ contains val = amg_poly_lottes_beta_ case('POLY_NEW') val = amg_poly_new_ + case('POLY_RHO_EST_POWER') + val = amg_poly_rho_est_power_ case('A_NORMI') val = amg_max_norm_ case('USER_CHOICE') diff --git a/amgprec/amg_d_poly_smoother.f90 b/amgprec/amg_d_poly_smoother.f90 index 7e444550..2d0ac1e1 100644 --- a/amgprec/amg_d_poly_smoother.f90 +++ b/amgprec/amg_d_poly_smoother.f90 @@ -59,9 +59,11 @@ module amg_d_poly_smoother ! class(amg_d_base_solver_type), allocatable :: sv ! integer(psb_ipk_) :: pdegree, variant + integer(psb_ipk_) :: rho_estimate=amg_poly_rho_est_power_ + integer(psb_ipk_) :: rho_estimate_iterations=10 type(psb_dspmat_type), pointer :: pa => null() real(psb_dpk_), allocatable :: poly_beta(:) - real(psb_dpk_) :: rho_ba + real(psb_dpk_) :: rho_ba = -done contains procedure, pass(sm) :: apply_v => amg_d_poly_smoother_apply_vect !!$ procedure, pass(sm) :: apply_a => amg_d_poly_smoother_apply @@ -317,10 +319,11 @@ contains ! ! Default: BJAC with no residual check ! - sm%pdegree = 1 - sm%rho_ba = dzero - sm%variant = amg_poly_lottes_ - + sm%pdegree = 1 + sm%rho_ba = -done + sm%variant = amg_poly_lottes_ + sm%rho_estimate = amg_poly_rho_est_power_ + sm%rho_estimate_iterations = 20 if (allocated(sm%sv)) then call sm%sv%default() end if diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 index 2158e3ae..d668bac0 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 @@ -98,44 +98,50 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) goto 9999 end if - if (.false.) then - 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 +!!$ if (.false.) then +!!$ 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 + if (sm%rho_ba <= dzero) then + select case(sm%rho_estimate) + case(amg_poly_rho_est_power_) + block + 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,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,sm%rho_estimate_iterations + 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 + end block + case default + write(0,*) ' Unknown algorithm for RHO(BA) estimate, defaulting to a value of 1.0 ' + sm%rho_ba = done end select - else - block - 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,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 if (debug_level >= psb_debug_outer_) & diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 index 0daa387b..e61b09e3 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 @@ -55,6 +55,8 @@ subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx) 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_d_base_smoother_type%set(what,val,info,idx=idx) end select diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 index 0b116deb..916fb5e6 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 @@ -64,6 +64,21 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx) 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_d_base_smoother_type%set(what,val,info,idx=idx) end select diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 index de308a8e..f1987a7b 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 @@ -54,8 +54,13 @@ subroutine amg_d_poly_smoother_csetr(sm,what,val,info,idx) call psb_erractionsave(err_act) select case(psb_toupper(what)) - case('RHO_BA') - sm%rho_ba = val + case('POLY_RHO_BA') + if ((dzero