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

Added timers and removed unuseful spmm

Browse Source
l1aggregation
Cirdans-Home 8 months ago
parent e83bde6896
commit 74dccb6c44
10 changed files with 127 additions and 203 deletions
Show Stats Download Patch File Download Diff File

6
amgprec/amg_d_matchboxp_mod.f90
Unescape Escape View File

@ -145,7 +145,7 @@ contains
logical, parameter :: dump=.false., debug=.false., dump_mate=.false., & logical, parameter :: dump=.false., debug=.false., dump_mate=.false., &
& debug_ilaggr=.false., debug_sync=.false., debug_mate=.false. & debug_ilaggr=.false., debug_sync=.false., debug_mate=.false.
integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1 integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1
logical, parameter :: do_timings=.false. logical, parameter :: do_timings=.true.
integer, parameter :: ilaggr_neginit=-1, ilaggr_nonlocal=-2 integer, parameter :: ilaggr_neginit=-1, ilaggr_nonlocal=-2
ictxt = desc_a%get_ctxt() ictxt = desc_a%get_ctxt()
@ -608,7 +608,7 @@ contains
logical, parameter :: old_style=.false., sort_minp=.true. logical, parameter :: old_style=.false., sort_minp=.true.
character(len=40) :: name='build_matching', fname character(len=40) :: name='build_matching', fname
integer(psb_ipk_), save :: idx_cmboxp=-1, idx_bldahat=-1, idx_phase2=-1, idx_phase3=-1 integer(psb_ipk_), save :: idx_cmboxp=-1, idx_bldahat=-1, idx_phase2=-1, idx_phase3=-1
logical, parameter :: do_timings=.false. logical, parameter :: do_timings=.true.
ictxt = desc_a%get_ctxt() ictxt = desc_a%get_ctxt()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -810,7 +810,7 @@ contains
character(len=80) :: aname character(len=80) :: aname
real(psb_dpk_), parameter :: eps=epsilon(1.d0) real(psb_dpk_), parameter :: eps=epsilon(1.d0)
integer(psb_ipk_), save :: idx_glbt=-1, idx_phase1=-1, idx_phase2=-1 integer(psb_ipk_), save :: idx_glbt=-1, idx_phase1=-1, idx_phase2=-1
logical, parameter :: do_timings=.false. logical, parameter :: do_timings=.true.
logical, parameter :: debug_symmetry = .false., check_size=.false. logical, parameter :: debug_symmetry = .false., check_size=.false.
logical, parameter :: unroll_logtrans=.false. logical, parameter :: unroll_logtrans=.false.

2
amgprec/impl/aggregator/amg_d_parmatch_aggregator_tprol.F90
Unescape Escape View File

@ -88,7 +88,7 @@ subroutine amg_d_parmatch_aggregator_build_tprol(ag,parms,ag_data,&
type(psb_ldspmat_type) :: tmp_prol, tmp_pg, tmp_restr type(psb_ldspmat_type) :: tmp_prol, tmp_pg, tmp_restr
type(psb_desc_type) :: tmp_desc_ac, tmp_desc_ax, tmp_desc_p type(psb_desc_type) :: tmp_desc_ac, tmp_desc_ax, tmp_desc_p
integer(psb_ipk_), save :: idx_mboxp=-1, idx_spmmbld=-1, idx_sweeps_mult=-1 integer(psb_ipk_), save :: idx_mboxp=-1, idx_spmmbld=-1, idx_sweeps_mult=-1
logical, parameter :: dump=.false., do_timings=.false., debug=.false., & logical, parameter :: dump=.false., do_timings=.true., debug=.false., &
& dump_prol_restr=.false. & dump_prol_restr=.false.
name='d_parmatch_tprol' name='d_parmatch_tprol'

2
amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_inner.F90
Unescape Escape View File

@ -131,7 +131,7 @@ subroutine amg_d_parmatch_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,&
& nzt, naggrm1, naggrp1, i, k & nzt, naggrm1, naggrp1, i, k
integer(psb_lpk_), allocatable :: ia(:),ja(:) integer(psb_lpk_), allocatable :: ia(:),ja(:)
!integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza, nrpsave, ncpsave, nzpsave !integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza, nrpsave, ncpsave, nzpsave
logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false.
integer(psb_ipk_), save :: idx_spspmm=-1, idx_prolcnv=-1, idx_proltrans=-1, idx_asb=-1 integer(psb_ipk_), save :: idx_spspmm=-1, idx_prolcnv=-1, idx_proltrans=-1, idx_asb=-1
name='amg_parmatch_spmm_bld_inner' name='amg_parmatch_spmm_bld_inner'

2
amgprec/impl/aggregator/amg_d_ptap_bld.f90
Unescape Escape View File

@ -486,7 +486,7 @@ subroutine amg_ld_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, &
& nzt, naggrm1, naggrp1, i, k & nzt, naggrm1, naggrp1, i, k
integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza
logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false.
integer(psb_ipk_), save :: idx_spspmm=-1 integer(psb_ipk_), save :: idx_spspmm=-1
name='amg_ptap_bld' name='amg_ptap_bld'

2
amgprec/impl/aggregator/amg_d_soc2_map_bld.F90
Unescape Escape View File

@ -104,7 +104,7 @@ subroutine amg_d_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
integer(psb_ipk_), save :: idx_soc2_p1=-1, idx_soc2_p2=-1, idx_soc2_p3=-1 integer(psb_ipk_), save :: idx_soc2_p1=-1, idx_soc2_p2=-1, idx_soc2_p3=-1
integer(psb_ipk_), save :: idx_soc2_p0=-1 integer(psb_ipk_), save :: idx_soc2_p0=-1
logical, parameter :: do_timings=.false. logical, parameter :: do_timings=.true.
info=psb_success_ info=psb_success_
name = 'amg_soc2_map_bld' name = 'amg_soc2_map_bld'

34
amgprec/impl/amg_dmlprec_aply.f90
Unescape Escape View File

@ -591,8 +591,6 @@ contains
integer(psb_ipk_) :: nlev, ilev, sweeps integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post logical :: pre, post
character(len=20) :: name character(len=20) :: name
logical, parameter :: do_timings=.true.
integer(psb_ipk_), save :: ml_mlt_smth=-1, ml_mlt_rp=-1, ml_mlt_rsd=-1
name = 'inner_inner_mult' name = 'inner_inner_mult'
info = psb_success_ info = psb_success_
@ -610,12 +608,6 @@ contains
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level write(debug_unit,*) me,' inner_mult at level ',level
end if end if
if ((do_timings).and.(ml_mlt_smth==-1)) &
& ml_mlt_smth = psb_get_timer_idx("ML-MLT: smoother ")
if ((do_timings).and.(ml_mlt_rp==-1)) &
& ml_mlt_rp = psb_get_timer_idx("ML-MLT: RestProl")
if ((do_timings).and.(ml_mlt_rsd==-1)) &
& ml_mlt_rsd = psb_get_timer_idx("ML-MLT: Residual")
sweeps_post = p%precv(level)%parms%sweeps_post sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre sweeps_pre = p%precv(level)%parms%sweeps_pre
@ -631,7 +623,7 @@ contains
! Apply the first smoother ! Apply the first smoother
! The residual has been prepared before the recursive call. ! The residual has been prepared before the recursive call.
! !
if (do_timings) call psb_tic(ml_mlt_smth)
if (pre) then if (pre) then
if (me >=0) then if (me >=0) then
!!$ write(0,*) me,'Applying smoother pre ', level !!$ write(0,*) me,'Applying smoother pre ', level
@ -654,13 +646,10 @@ contains
end if end if
end if end if
endif endif
if (do_timings) call psb_toc(ml_mlt_smth)
! !
! Compute the residual for next level and call recursively ! Compute the residual for next level and call recursively
! !
if (pre) then if (pre) then
if (do_timings) call psb_tic(ml_mlt_rsd)
call psb_geaxpby(done,vx2l,& call psb_geaxpby(done,vx2l,&
& dzero,vty,& & dzero,vty,&
& base_desc,info) & base_desc,info)
@ -668,9 +657,6 @@ contains
if (info == psb_success_) call psb_spmm(-done,base_a,& if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,& & vy2l,done,vty,&
& base_desc,info,work=work,trans=trans) & base_desc,info,work=work,trans=trans)
if (do_timings) call psb_toc(ml_mlt_rsd)
if (do_timings) call psb_tic(ml_mlt_rp)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue') & a_err='Error during residue')
@ -685,9 +671,7 @@ contains
& a_err='Error during restriction') & a_err='Error during restriction')
goto 9999 goto 9999
end if end if
if (do_timings) call psb_toc(ml_mlt_rp)
else else
if (do_timings) call psb_tic(ml_mlt_rp)
! Shortcut: just transfer x2l. ! Shortcut: just transfer x2l.
call p%precv(level+1)%map_rstr(done,vx2l,& call p%precv(level+1)%map_rstr(done,vx2l,&
& dzero,p%precv(level+1)%wrk%vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,&
@ -698,7 +682,6 @@ contains
& a_err='Error during restriction') & a_err='Error during restriction')
goto 9999 goto 9999
end if end if
if (do_timings) call psb_toc(ml_mlt_rp)
endif endif
call inner_ml_aply(level+1,p,trans,work,info) call inner_ml_aply(level+1,p,trans,work,info)
@ -706,12 +689,10 @@ contains
! !
! Apply the prolongator ! Apply the prolongator
! !
if (do_timings) call psb_tic(ml_mlt_rp)
call p%precv(level+1)%map_prol(done,& call p%precv(level+1)%map_prol(done,&
& p%precv(level+1)%wrk%vy2l,done,vy2l,& & p%precv(level+1)%wrk%vy2l,done,vy2l,&
& info,work=work,& & info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
if (do_timings) call psb_toc(ml_mlt_rp)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -719,7 +700,7 @@ contains
end if end if
if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then if (p%precv(level)%parms%ml_cycle == amg_wcycle_ml_) then
if (do_timings) call psb_tic(ml_mlt_rsd)
if (me >=0) then if (me >=0) then
call psb_geaxpby(done,vx2l, dzero,vty,& call psb_geaxpby(done,vx2l, dzero,vty,&
& base_desc,info) & base_desc,info)
@ -727,13 +708,10 @@ contains
& vy2l,done,vty,& & vy2l,done,vty,&
& base_desc,info,work=work,trans=trans) & base_desc,info,work=work,trans=trans)
end if end if
if (do_timings) call psb_toc(ml_mlt_rsd)
if (do_timings) call psb_tic(ml_mlt_rp)
if (info == psb_success_) & if (info == psb_success_) &
& call p%precv(level+1)%map_rstr(done,vty,& & call p%precv(level+1)%map_rstr(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,info,work=work,& & dzero,p%precv(level+1)%wrk%vx2l,info,work=work,&
& vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1))
if (do_timings) call psb_toc(ml_mlt_rp)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction') & a_err='Error during W-cycle restriction')
@ -742,12 +720,10 @@ contains
call inner_ml_aply(level+1,p,trans,work,info) call inner_ml_aply(level+1,p,trans,work,info)
if (do_timings) call psb_tic(ml_mlt_rp)
if (info == psb_success_) call p%precv(level+1)%map_prol(done, & if (info == psb_success_) call p%precv(level+1)%map_prol(done, &
& p%precv(level+1)%wrk%vy2l,done,vy2l,& & p%precv(level+1)%wrk%vy2l,done,vy2l,&
& info,work=work,& & info,work=work,&
& vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1))
if (do_timings) call psb_toc(ml_mlt_rp)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -760,7 +736,6 @@ contains
if (post) then if (post) then
if (me >=0) then if (me >=0) then
if (do_timings) call psb_tic(ml_mlt_rsd)
call psb_geaxpby(done,vx2l,& call psb_geaxpby(done,vx2l,&
& dzero,vty,& & dzero,vty,&
& base_desc,info) & base_desc,info)
@ -772,9 +747,7 @@ contains
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
if (do_timings) call psb_toc(ml_mlt_rsd)
if (do_timings) call psb_tic(ml_mlt_smth)
! !
! Apply the second smoother ! Apply the second smoother
! !
@ -789,7 +762,6 @@ contains
& vty,done,vy2l, base_desc, trans,& & vty,done,vy2l, base_desc, trans,&
& sweeps,work,wv,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (do_timings) call psb_toc(ml_mlt_smth)
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -802,14 +774,12 @@ contains
else if (level == nlev) then else if (level == nlev) then
!!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal() !!$ write(0,*) me,'Applying smoother at top level ',psb_errstatus_fatal()
if (do_timings) call psb_tic(ml_mlt_smth)
if (me >=0) then if (me >=0) then
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,base_desc, trans,& & vx2l,dzero,vy2l,base_desc, trans,&
& sweeps,work,wv,info) & sweeps,work,wv,info)
end if end if
if (do_timings) call psb_toc(ml_mlt_smth)
!!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal() !!$ write(0,*) me,' Done applying smoother at top level ',psb_errstatus_fatal()
else else

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

@ -55,6 +55,10 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_d_vect_type),intent(inout), optional :: initu type(psb_d_vect_type),intent(inout), optional :: initu
! Timers
logical, parameter :: do_timings=.true.
integer(psb_ipk_), save :: poly_1=-1, poly_2=-1, poly_3=-1
integer(psb_ipk_), save :: poly_mv=-1, poly_sv=-1, poly_vect=-1
! !
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_d_vect_type) :: tx, ty, tz, r type(psb_d_vect_type) :: tx, ty, tz, r
@ -64,9 +68,6 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character :: trans_, init_ character :: trans_, init_
real(psb_dpk_) :: res, resdenum real(psb_dpk_) :: res, resdenum
character(len=20) :: name='d_poly_smoother_apply_v' character(len=20) :: name='d_poly_smoother_apply_v'
logical, parameter :: do_timings=.true.
integer(psb_ipk_), save :: poly_1=-1, poly_2=-1, poly_3=-1
integer(psb_ipk_), save :: poly_mv=-1, poly_sv=-1, poly_vect=-1
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -146,35 +147,33 @@ 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, sm%pdegree do i=1, sm%pdegree-1
! 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') if (do_timings) call psb_tic(poly_sv)
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
if (do_timings) call psb_toc(poly_sv)
cz = (2*i*done-3)/(2*i*done+done) cz = (2*i*done-3)/(2*i*done+done)
cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
if (.false.) then if (do_timings) call psb_tic(poly_vect)
! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1
call psb_geaxpby(cr,ty,cz,tz,desc_data,info) if (do_timings) call psb_toc(poly_vect)
! r_k = b-Ax_k = x -A tx if (do_timings) call psb_tic(poly_mv)
call psb_geaxpby(done,tz,done,tx,desc_data,info) call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
else if (do_timings) call psb_toc(poly_mv)
call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
end if
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 LOTTES',i,res
! x_k = x_{k-1} + z_k
end do end do
if (do_timings) call psb_tic(poly_sv)
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
if (do_timings) call psb_toc(poly_sv)
cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
if (do_timings) call psb_tic(poly_vect)
call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
if (do_timings) call psb_toc(poly_vect)
end block end block
if (do_timings) call psb_toc(poly_1) if (do_timings) call psb_toc(poly_1)
case(amg_poly_lottes_beta_) case(amg_poly_lottes_beta_)
if (do_timings) call psb_tic(poly_2) if (do_timings) call psb_tic(poly_2)
block block
real(psb_dpk_) :: cz, cr real(psb_dpk_) :: cz, cr
! b == x ! b == x
@ -188,42 +187,35 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
end if end if
do i=1, sm%pdegree do i=1, sm%pdegree-1
! B r_{k-1} ! B r_{k-1}
if (do_timings) call psb_tic(poly_sv)
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')
if (do_timings) call psb_toc(poly_sv)
cz = (2*i*done-3)/(2*i*done+done) cz = (2*i*done-3)/(2*i*done+done)
cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
if (.false.) then if (do_timings) call psb_tic(poly_vect)
! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} call psb_abgdxyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
call psb_geaxpby(cr,ty,cz,tz,desc_data,info) if (do_timings) call psb_toc(poly_vect)
! r_k = b-Ax_k = x -A tx if (do_timings) call psb_tic(poly_mv)
call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info) call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
else if (do_timings) call psb_toc(poly_mv)
call psb_abgdxyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
end if
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 LOTTES_BETA',i,res
! x_k = x_{k-1} + z_k
end do end do
call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
if (do_timings) call psb_tic(poly_vect)
call psb_abgdxyz(cr,cz,sm%poly_beta(sm%pdegree),done,ty,tz,tx,desc_data,info)
if (do_timings) call psb_toc(poly_vect)
end block end block
if (do_timings) call psb_toc(poly_2) if (do_timings) call psb_toc(poly_2)
case(amg_poly_new_) case(amg_poly_new_)
if (do_timings) call psb_tic(poly_3) if (do_timings) call psb_tic(poly_3)
block block
real(psb_dpk_) :: sigma, theta, delta, rho_old, rho real(psb_dpk_) :: sigma, theta, delta, rho_old, rho
! b == x ! b == x
! x == tx ! x == tx
! !
sm%rho_ba = 1.12d0
!write(0,*) 'Parameter: ',sm%cf_a,sm%rho_ba
theta = (done+sm%cf_a)/2 theta = (done+sm%cf_a)/2
delta = (done-sm%cf_a)/2 delta = (done-sm%cf_a)/2
@ -232,21 +224,15 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (do_timings) call psb_tic(poly_sv) if (do_timings) call psb_tic(poly_sv)
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')
if (do_timings) call psb_toc(poly_sv) if (do_timings) call psb_toc(poly_sv)
if (do_timings) call psb_tic(poly_vect)
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)
!write(0,*) 'POLY_NEW Iteration',0,' :',psb_genrm2(r,desc_data,info) if (do_timings) call psb_tic(poly_vect)
if (.false.) then call psb_abgdxyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info)
call psb_geaxpby(done,tz,done,tx,desc_data,info)
else
call psb_abgdxyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
end if
if (do_timings) call psb_toc(poly_vect) if (do_timings) call psb_toc(poly_vect)
! tz == d ! tz == d
do i=1, sm%pdegree-1 do i=1, sm%pdegree-1
! !
! !
! r_{k-1} = r_k - (1/rho(BA)) B A d_k ! r_{k-1} = r_k - (1/rho(BA)) B A d_k
if (do_timings) call psb_tic(poly_mv) if (do_timings) call psb_tic(poly_mv)
call psb_spmm(done,sm%pa,tz,dzero,ty,desc_data,info,work=aux,trans=trans_) call psb_spmm(done,sm%pa,tz,dzero,ty,desc_data,info,work=aux,trans=trans_)
@ -254,54 +240,16 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (do_timings) call psb_tic(poly_sv) if (do_timings) call psb_tic(poly_sv)
call sm%sv%apply(-(done/sm%rho_ba),ty,done,r,desc_data,trans_,aux,wv(5:),info,init='Z') call sm%sv%apply(-(done/sm%rho_ba),ty,done,r,desc_data,trans_,aux,wv(5:),info,init='Z')
if (do_timings) call psb_toc(poly_sv) if (do_timings) call psb_toc(poly_sv)
if (do_timings) call psb_tic(poly_vect)
!write(0,*) 'POLY_NEW Iteration',i,' :',psb_genrm2(r,desc_data,info)
! !
! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1} ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
rho = done/(2*sigma - rho_old) rho = done/(2*sigma - rho_old)
if (.false.) then if (do_timings) call psb_tic(poly_vect)
call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info) call psb_abgdxyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
call psb_geaxpby(done,tz,done,tx,desc_data,info)
else
call psb_abgdxyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
end if
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Polynomial smoother NEW ',i,res
! x_k = x_{k-1} + z_k
rho_old = rho
if (do_timings) call psb_toc(poly_vect) if (do_timings) call psb_toc(poly_vect)
end do
end block
if (do_timings) call psb_toc(poly_3)
case(amg_poly_dbg_)
block
real(psb_dpk_) :: sigma, theta, delta, rho_old, rho
! b == x
! x == tx
!
write(0,*) 'Parameter: ',sm%cf_a
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,ty,dzero,r,desc_data,info)
call psb_geaxpby(done/theta,r,dzero,tz,desc_data,info)
write(0,*) 'POLY_DBG Iteration',0,' :',psb_genrm2(r,desc_data,info)
do i=1, sm%pdegree
call psb_geaxpby(done,tz,done,tx,desc_data,info)
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')
write(0,*) 'POLY_DBG Iteration',i,' :',psb_genrm2(r,desc_data,info)
rho = done/(2*sigma - rho_old)
call psb_geaxpby((2*rho/delta),r,rho*rho_old,tz,desc_data,info)
rho_old = rho rho_old = rho
end do end do
end block end block
if (do_timings) call psb_toc(poly_3)
case default case default
info=psb_err_internal_error_ info=psb_err_internal_error_
call psb_errpush(info,name,& call psb_errpush(info,name,&

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

@ -87,14 +87,10 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='invalid sm%degree for poly_beta') & a_err='invalid sm%degree for poly_beta')
goto 9999 goto 9999
end if end if
case(amg_poly_new_, amg_poly_dbg_) case(amg_poly_new_)
if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
!Ok !Ok
!!$ write(0,*) 'Vector: '
!!$ do i=1,size(amg_d_poly_a_vect)
!!$ write(0,*) i,amg_d_poly_a_vect(i)
!!$ end do
sm%cf_a = amg_d_poly_a_vect(sm%pdegree) sm%cf_a = amg_d_poly_a_vect(sm%pdegree)
else else
info = psb_err_internal_error_ info = psb_err_internal_error_

2
amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90
Unescape Escape View File

@ -58,7 +58,7 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx)
sm%pdegree = val sm%pdegree = val
case('POLY_VARIANT') case('POLY_VARIANT')
select case(val) select case(val)
case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_,amg_poly_dbg_) case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_)
sm%variant = val sm%variant = val
case default case default
write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val

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

@ -55,6 +55,10 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu type(psb_s_vect_type),intent(inout), optional :: initu
! Timers
logical, parameter :: do_timings=.true.
integer(psb_ipk_), save :: poly_1=-1, poly_2=-1, poly_3=-1
integer(psb_ipk_), save :: poly_mv=-1, poly_sv=-1, poly_vect=-1
! !
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_s_vect_type) :: tx, ty, tz, r type(psb_s_vect_type) :: tx, ty, tz, r
@ -93,6 +97,18 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
if ((do_timings).and.(poly_1==-1)) &
& poly_1 = psb_get_timer_idx("POLY: Chebychev4")
if ((do_timings).and.(poly_2==-1)) &
& poly_2 = psb_get_timer_idx("POLY: OptChebychev4")
if ((do_timings).and.(poly_3==-1)) &
& poly_3 = psb_get_timer_idx("POLY: OptChebychev1")
if ((do_timings).and.(poly_mv==-1)) &
& poly_mv = psb_get_timer_idx("POLY: spMV")
if ((do_timings).and.(poly_vect==-1)) &
& poly_vect = psb_get_timer_idx("POLY: Vectors")
if ((do_timings).and.(poly_sv==-1)) &
& poly_sv = psb_get_timer_idx("POLY: solver")
n_row = desc_data%get_local_rows() n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols() n_col = desc_data%get_local_cols()
@ -125,38 +141,39 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
select case(sm%variant) select case(sm%variant)
case(amg_poly_lottes_) case(amg_poly_lottes_)
if (do_timings) call psb_tic(poly_1)
block block
real(psb_spk_) :: cz, cr real(psb_spk_) :: cz, cr
! b == x ! b == x
! x == tx ! x == tx
! !
do i=1, sm%pdegree do i=1, sm%pdegree-1
! 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') if (do_timings) call psb_tic(poly_sv)
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
if (do_timings) call psb_toc(poly_sv)
cz = (2*i*sone-3)/(2*i*sone+sone) cz = (2*i*sone-3)/(2*i*sone+sone)
cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba) cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
if (.false.) then if (do_timings) call psb_tic(poly_vect)
! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1
call psb_geaxpby(cr,ty,cz,tz,desc_data,info) if (do_timings) call psb_toc(poly_vect)
! r_k = b-Ax_k = x -A tx if (do_timings) call psb_tic(poly_mv)
call psb_geaxpby(sone,tz,sone,tx,desc_data,info) call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
else if (do_timings) call psb_toc(poly_mv)
call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
end if
if (.false.) then
call psb_geaxpby(sone,x,szero,r,desc_data,info)
call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_)
else
call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
end if
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Polynomial smoother LOTTES',i,res
! x_k = x_{k-1} + z_k
end do end do
if (do_timings) call psb_tic(poly_sv)
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
if (do_timings) call psb_toc(poly_sv)
cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
if (do_timings) call psb_tic(poly_vect)
call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
if (do_timings) call psb_toc(poly_vect)
end block end block
if (do_timings) call psb_toc(poly_1)
case(amg_poly_lottes_beta_) case(amg_poly_lottes_beta_)
if (do_timings) call psb_tic(poly_2)
block block
real(psb_spk_) :: cz, cr real(psb_spk_) :: cz, cr
! b == x ! b == x
@ -170,32 +187,30 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
end if end if
do i=1, sm%pdegree do i=1, sm%pdegree-1
! B r_{k-1} ! B r_{k-1}
if (do_timings) call psb_tic(poly_sv)
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')
if (do_timings) call psb_toc(poly_sv)
cz = (2*i*sone-3)/(2*i*sone+sone) cz = (2*i*sone-3)/(2*i*sone+sone)
cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba) cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
if (.false.) then if (do_timings) call psb_tic(poly_vect)
! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} call psb_abgdxyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
call psb_geaxpby(cr,ty,cz,tz,desc_data,info) if (do_timings) call psb_toc(poly_vect)
! r_k = b-Ax_k = x -A tx if (do_timings) call psb_tic(poly_mv)
call psb_geaxpby(sm%poly_beta(i),tz,sone,tx,desc_data,info) call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
else if (do_timings) call psb_toc(poly_mv)
call psb_abgdxyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
end if
if (.false.) then
call psb_geaxpby(sone,x,szero,r,desc_data,info)
call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_)
else
call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
end if
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Polynomial smoother LOTTES_BETA',i,res
! x_k = x_{k-1} + z_k
end do end do
call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
if (do_timings) call psb_tic(poly_vect)
call psb_abgdxyz(cr,cz,sm%poly_beta(sm%pdegree),sone,ty,tz,tx,desc_data,info)
if (do_timings) call psb_toc(poly_vect)
end block end block
if (do_timings) call psb_toc(poly_2)
case(amg_poly_new_) case(amg_poly_new_)
if (do_timings) call psb_tic(poly_3)
block block
real(psb_spk_) :: sigma, theta, delta, rho_old, rho real(psb_spk_) :: sigma, theta, delta, rho_old, rho
! b == x ! b == x
@ -206,40 +221,35 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
delta = (sone-sm%cf_a)/2 delta = (sone-sm%cf_a)/2
sigma = theta/delta sigma = theta/delta
rho_old = sone/sigma rho_old = sone/sigma
if (do_timings) call psb_tic(poly_sv)
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')
if (do_timings) call psb_toc(poly_sv)
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)
if (.false.) then if (do_timings) call psb_tic(poly_vect)
call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info) call psb_abgdxyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
call psb_geaxpby(sone,tz,sone,tx,desc_data,info) if (do_timings) call psb_toc(poly_vect)
else
call psb_abgdxyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
end if
! tz == d ! tz == d
do i=1, sm%pdegree-1 do i=1, sm%pdegree-1
! !
! !
! r_{k-1} = r_k - (1/rho(BA)) B A d_k ! r_{k-1} = r_k - (1/rho(BA)) B A d_k
if (do_timings) call psb_tic(poly_mv)
call psb_spmm(sone,sm%pa,tz,szero,ty,desc_data,info,work=aux,trans=trans_) call psb_spmm(sone,sm%pa,tz,szero,ty,desc_data,info,work=aux,trans=trans_)
if (do_timings) call psb_toc(poly_mv)
if (do_timings) call psb_tic(poly_sv)
call sm%sv%apply(-(sone/sm%rho_ba),ty,sone,r,desc_data,trans_,aux,wv(5:),info,init='Z') call sm%sv%apply(-(sone/sm%rho_ba),ty,sone,r,desc_data,trans_,aux,wv(5:),info,init='Z')
if (do_timings) call psb_toc(poly_sv)
! !
! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1} ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
rho = sone/(2*sigma - rho_old) rho = sone/(2*sigma - rho_old)
if (.false.) then if (do_timings) call psb_tic(poly_vect)
call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info) call psb_abgdxyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
call psb_geaxpby(sone,tz,sone,tx,desc_data,info) if (do_timings) call psb_toc(poly_vect)
else
call psb_abgdxyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
end if
!!$ res = psb_genrm2(r,desc_data,info)
!!$ write(0,*) 'Polynomial smoother NEW ',i,res
! x_k = x_{k-1} + z_k
rho_old = rho rho_old = rho
end do end do
end block end block
if (do_timings) call psb_toc(poly_3)
case default case default
info=psb_err_internal_error_ info=psb_err_internal_error_
call psb_errpush(info,name,& call psb_errpush(info,name,&

Write Preview
Loading…
Cancel
Save