|
|
@ -36,7 +36,7 @@
|
|
|
|
!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
|
|
|
subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
|
|
|
& sweeps,work,wv,info,init,initu)
|
|
|
|
& sweeps,work,wv,info,init,initu)
|
|
|
|
|
|
|
|
|
|
|
|
use psb_base_mod
|
|
|
|
use psb_base_mod
|
|
|
|
use amg_d_diag_solver
|
|
|
|
use amg_d_diag_solver
|
|
|
@ -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
|
|
|
@ -92,7 +96,19 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
|
|
|
call psb_errpush(info,name)
|
|
|
|
call psb_errpush(info,name)
|
|
|
|
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_d_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_dpk_) :: cz, cr
|
|
|
|
real(psb_dpk_) :: 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(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)
|
|
|
|
|
|
|
|
|
|
|
|
case(amg_poly_lottes_beta_)
|
|
|
|
case(amg_poly_lottes_beta_)
|
|
|
|
|
|
|
|
if (do_timings) call psb_tic(poly_2)
|
|
|
|
block
|
|
|
|
block
|
|
|
|
real(psb_dpk_) :: cz, cr
|
|
|
|
real(psb_dpk_) :: cz, cr
|
|
|
|
! b == x
|
|
|
|
! b == x
|
|
|
@ -170,32 +187,30 @@ 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)
|
|
|
|
case(amg_poly_new_)
|
|
|
|
case(amg_poly_new_)
|
|
|
|
|
|
|
|
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
|
|
|
@ -206,40 +221,35 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
|
|
|
|
delta = (done-sm%cf_a)/2
|
|
|
|
delta = (done-sm%cf_a)/2
|
|
|
|
sigma = theta/delta
|
|
|
|
sigma = theta/delta
|
|
|
|
rho_old = done/sigma
|
|
|
|
rho_old = done/sigma
|
|
|
|
|
|
|
|
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)
|
|
|
|
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)
|
|
|
|
if (.false.) then
|
|
|
|
if (do_timings) call psb_tic(poly_vect)
|
|
|
|
call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info)
|
|
|
|
call psb_abgdxyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
|
|
|
|
call psb_geaxpby(done,tz,done,tx,desc_data,info)
|
|
|
|
if (do_timings) call psb_toc(poly_vect)
|
|
|
|
else
|
|
|
|
|
|
|
|
call psb_abgdxyz((done/theta),dzero,done,done,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(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_)
|
|
|
|
|
|
|
|
if (do_timings) call psb_toc(poly_mv)
|
|
|
|
|
|
|
|
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)
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! 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)
|
|
|
|
if (do_timings) call psb_toc(poly_vect)
|
|
|
|
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
|
|
|
|
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,&
|
|
|
|