fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed build and apply to actually use degree

Browse Source
Poly-novrl
Cirdans-Home 12 months ago
parent 54d608d2dd
commit ea8974f88c
6 changed files with 100 additions and 102 deletions
Show Stats Download Patch File Download Diff File

18
amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
type(psb_d_vect_type),intent(inout) :: y type(psb_d_vect_type),intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps ! this is ignored here, the polynomial degree dictates the value
real(psb_dpk_),target, intent(inout) :: work(:) real(psb_dpk_),target, intent(inout) :: work(:)
type(psb_d_vect_type),intent(inout) :: wv(:) type(psb_d_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -115,7 +115,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& a_err='invalid wv size in smoother_apply') & a_err='invalid wv size in smoother_apply')
goto 9999 goto 9999
end if end if
sm%pdegree = sweeps
associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4)) associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4))
call psb_geaxpby(done,x,dzero,r,desc_data,info) call psb_geaxpby(done,x,dzero,r,desc_data,info)
@ -130,7 +130,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! b == x ! b == x
! x == tx ! x == tx
! !
do i=1, sweeps do i=1, sm%pdegree
! B r_{k-1} ! B r_{k-1}
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') 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) cz = (2*i*done-3)/(2*i*done+done)
@ -159,14 +159,14 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! x == tx ! x == tx
! !
if (allocated(sm%poly_beta)) then if (allocated(sm%poly_beta)) then
if (size(sm%poly_beta) /= sweeps) deallocate(sm%poly_beta) if (size(sm%poly_beta) /= sm%pdegree) deallocate(sm%poly_beta)
end if end if
if (.not.allocated(sm%poly_beta)) then if (.not.allocated(sm%poly_beta)) then
call psb_realloc(sweeps,sm%poly_beta,info) call psb_realloc(sm%pdegree,sm%poly_beta,info)
sm%poly_beta(1:sweeps) = amg_d_poly_beta_mat(1:sweeps,sweeps) sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
end if end if
do i=1, sweeps do i=1, sm%pdegree
! B r_{k-1} ! B r_{k-1}
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') 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) cz = (2*i*done-3)/(2*i*done+done)
@ -193,7 +193,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! b == x ! b == x
! x == tx ! x == tx
! !
sm%cf_a = amg_d_poly_a_vect(sweeps)
theta = (done+sm%cf_a)/2 theta = (done+sm%cf_a)/2
delta = (done-sm%cf_a)/2 delta = (done-sm%cf_a)/2
@ -203,7 +203,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/sm%rho_ba),ty,dzero,r,desc_data,info)
call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info) call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info)
! tz == d ! tz == d
do i=1, sweeps do i=1, sm%pdegree
! x_{k+1} = x_k + d_k ! x_{k+1} = x_k + d_k
call psb_geaxpby(done,tz,done,tx,desc_data,info) call psb_geaxpby(done,tz,done,tx,desc_data,info)
! !

7
amgprec/impl/smoother/amg_d_poly_smoother_bld.f90
Unescape Escape View File

@ -74,8 +74,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
n_col = desc_a%get_local_cols() n_col = desc_a%get_local_cols()
nrow_a = a%get_nrows() nrow_a = a%get_nrows()
nztota = a%get_nzeros() nztota = a%get_nzeros()
if (.false.) then select case(sm%variant)
select case(sm%variant)
case(amg_poly_lottes_) case(amg_poly_lottes_)
! do nothing ! do nothing
case(amg_poly_lottes_beta_) case(amg_poly_lottes_beta_)
@ -89,7 +88,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
goto 9999 goto 9999
end if end if
case(amg_poly_new_) case(amg_poly_new_)
write(psb_out_unit,*) "Nella fase di build sm%pdegree = ",sm%pdegree
if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
!Ok !Ok
sm%cf_a = amg_d_poly_a_vect(sm%pdegree) sm%cf_a = amg_d_poly_a_vect(sm%pdegree)
@ -106,7 +105,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='invalid sm%variant') & a_err='invalid sm%variant')
goto 9999 goto 9999
end select end select
end if
sm%pa => a sm%pa => a
if (.not.allocated(sm%sv)) then if (.not.allocated(sm%sv)) then
info = psb_err_internal_error_ info = psb_err_internal_error_

10
amgprec/impl/smoother/amg_d_poly_smoother_descr.f90
Unescape Escape View File

@ -79,18 +79,18 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix)
select case(sm%variant) select case(sm%variant)
case(amg_poly_lottes_) case(amg_poly_lottes_)
write(iout_,*) trim(prefix_), ' variant: ','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 write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
case(amg_poly_lottes_beta_) case(amg_poly_lottes_beta_)
write(iout_,*) trim(prefix_), ' variant: ','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 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_) case(amg_poly_new_)
write(iout_,*) trim(prefix_), ' variant: ','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_), ' rho_ba: ',sm%rho_ba
!write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a
case default case default
write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???' write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???'
end select end select

18
amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
type(psb_s_vect_type),intent(inout) :: y type(psb_s_vect_type),intent(inout) :: y
real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps ! this is ignored here, the polynomial degree dictates the value
real(psb_spk_),target, intent(inout) :: work(:) real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:) type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
@ -115,7 +115,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& a_err='invalid wv size in smoother_apply') & a_err='invalid wv size in smoother_apply')
goto 9999 goto 9999
end if end if
sm%pdegree = sweeps
associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4)) associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4))
call psb_geaxpby(sone,x,szero,r,desc_data,info) call psb_geaxpby(sone,x,szero,r,desc_data,info)
@ -130,7 +130,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! b == x ! b == x
! x == tx ! x == tx
! !
do i=1, sweeps do i=1, sm%pdegree
! B r_{k-1} ! B r_{k-1}
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') 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) cz = (2*i*sone-3)/(2*i*sone+sone)
@ -159,14 +159,14 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! x == tx ! x == tx
! !
if (allocated(sm%poly_beta)) then if (allocated(sm%poly_beta)) then
if (size(sm%poly_beta) /= sweeps) deallocate(sm%poly_beta) if (size(sm%poly_beta) /= sm%pdegree) deallocate(sm%poly_beta)
end if end if
if (.not.allocated(sm%poly_beta)) then if (.not.allocated(sm%poly_beta)) then
call psb_realloc(sweeps,sm%poly_beta,info) call psb_realloc(sm%pdegree,sm%poly_beta,info)
sm%poly_beta(1:sweeps) = amg_d_poly_beta_mat(1:sweeps,sweeps) sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
end if end if
do i=1, sweeps do i=1, sm%pdegree
! B r_{k-1} ! B r_{k-1}
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') 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) cz = (2*i*sone-3)/(2*i*sone+sone)
@ -193,7 +193,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! b == x ! b == x
! x == tx ! x == tx
! !
sm%cf_a = amg_d_poly_a_vect(sweeps)
theta = (sone+sm%cf_a)/2 theta = (sone+sm%cf_a)/2
delta = (sone-sm%cf_a)/2 delta = (sone-sm%cf_a)/2
@ -203,7 +203,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info) call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info)
call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info) call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info)
! tz == d ! tz == d
do i=1, sweeps do i=1, sm%pdegree
! x_{k+1} = x_k + d_k ! x_{k+1} = x_k + d_k
call psb_geaxpby(sone,tz,sone,tx,desc_data,info) call psb_geaxpby(sone,tz,sone,tx,desc_data,info)
! !

7
amgprec/impl/smoother/amg_s_poly_smoother_bld.f90
Unescape Escape View File

@ -74,8 +74,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
n_col = desc_a%get_local_cols() n_col = desc_a%get_local_cols()
nrow_a = a%get_nrows() nrow_a = a%get_nrows()
nztota = a%get_nzeros() nztota = a%get_nzeros()
if (.false.) then select case(sm%variant)
select case(sm%variant)
case(amg_poly_lottes_) case(amg_poly_lottes_)
! do nothing ! do nothing
case(amg_poly_lottes_beta_) case(amg_poly_lottes_beta_)
@ -89,7 +88,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
goto 9999 goto 9999
end if end if
case(amg_poly_new_) case(amg_poly_new_)
write(psb_out_unit,*) "Nella fase di build sm%pdegree = ",sm%pdegree
if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
!Ok !Ok
sm%cf_a = amg_d_poly_a_vect(sm%pdegree) sm%cf_a = amg_d_poly_a_vect(sm%pdegree)
@ -106,7 +105,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='invalid sm%variant') & a_err='invalid sm%variant')
goto 9999 goto 9999
end select end select
end if
sm%pa => a sm%pa => a
if (.not.allocated(sm%sv)) then if (.not.allocated(sm%sv)) then
info = psb_err_internal_error_ info = psb_err_internal_error_

10
amgprec/impl/smoother/amg_s_poly_smoother_descr.f90
Unescape Escape View File

@ -79,18 +79,18 @@ subroutine amg_s_poly_smoother_descr(sm,info,iout,coarse,prefix)
select case(sm%variant) select case(sm%variant)
case(amg_poly_lottes_) case(amg_poly_lottes_)
write(iout_,*) trim(prefix_), ' variant: ','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 write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba
case(amg_poly_lottes_beta_) case(amg_poly_lottes_beta_)
write(iout_,*) trim(prefix_), ' variant: ','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 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_) case(amg_poly_new_)
write(iout_,*) trim(prefix_), ' variant: ','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_), ' rho_ba: ',sm%rho_ba
!write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a
case default case default
write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???' write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???'
end select end select

Write Preview
Loading…
Cancel
Save