From 23aabd794d3031ae805dc2d36dd1a9c41862cd46 Mon Sep 17 00:00:00 2001 From: sfilippone Date: Mon, 20 Nov 2023 11:32:51 +0100 Subject: [PATCH] Defined new variant of polynomial smoother. --- amgprec/amg_d_poly_smoother.f90 | 2 +- .../amg_d_poly_smoother_apply_vect.f90 | 123 ++++++++++++------ .../impl/smoother/amg_d_poly_smoother_bld.f90 | 7 +- .../smoother/amg_d_poly_smoother_descr.f90 | 5 + 4 files changed, 90 insertions(+), 47 deletions(-) diff --git a/amgprec/amg_d_poly_smoother.f90 b/amgprec/amg_d_poly_smoother.f90 index 5ba83c24..a87fbb1b 100644 --- a/amgprec/amg_d_poly_smoother.f90 +++ b/amgprec/amg_d_poly_smoother.f90 @@ -63,7 +63,7 @@ module amg_d_poly_smoother integer(psb_ipk_) :: rho_estimate_iterations=10 type(psb_dspmat_type), pointer :: pa => null() real(psb_dpk_), allocatable :: poly_beta(:) - real(psb_dpk_), allocatable :: poly_a(:) + real(psb_dpk_) :: cf_a = dzero real(psb_dpk_) :: rho_ba = -done contains procedure, pass(sm) :: apply_v => amg_d_poly_smoother_apply_vect diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 index 044c4fb5..9b085f62 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 @@ -63,7 +63,6 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& integer(psb_ipk_) :: np, me, i, err_act character :: trans_, init_ real(psb_dpk_) :: res, resdenum - real(psb_dpk_) :: cz, cr character(len=20) :: name='d_poly_smoother_apply_v' call psb_erractionsave(err_act) @@ -131,62 +130,100 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& select case(sm%variant) case(amg_poly_lottes_) - ! 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(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 + block + real(psb_dpk_) :: cz, cr + ! 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(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 + ! x_k = x_{k-1} + z_k + end do + end block 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 + block + real(psb_dpk_) :: cz, cr + ! 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 - + ! x_k = x_{k-1} + z_k + end do + end block + case(amg_poly_new_) + block + real(psb_dpk_) :: sigma, theta, delta, rho_old, rho + ! b == x + ! x == tx + ! + theta = (done+sm%cf_a)/2 + delta = (done-sm%cf_a)/2 + sigma = theta/delta + rho_old = done/sigma + call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') + 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 + ! x_{k+1} = x_k + d_k + call psb_geaxpby(done,tz,done,tx,desc_data,info) + ! + ! r_{k-1} = r_k - (1/rho(BA)) B A d_k + call psb_spmm(done,sm%pa,tz,dzero,ty,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(-done,ty,done,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 = done/(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 diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 index 4c587f4e..a7f0a72c 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 @@ -90,15 +90,16 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) end if case(amg_poly_new_) - if ((1<=sm%pdegree).and.(sm%pdegree<=6)) then - call psb_realloc(sm%pdegree,sm%poly_a,info) - sm%poly_a(1:sm%pdegree) = amg_d_poly_a_vect(1:sm%pdegree) + 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,& diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 index 0607064d..97521259 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 @@ -88,6 +88,11 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) 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