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

mld2p4-2:

Browse Source 
mlprec/impl/mld_dmlprec_aply.f90
 mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90

Testing: make APPLY_VECT go through the iterative path, always; this also
handles the case SWEEPS=0
DMLPREC_APLY: activate Briggs-style option.
stopcriterion
Salvatore Filippone 9 years ago
parent 32e344aea4
commit e7492ad867
5 changed files with 192 additions and 155 deletions
Show Stats Download Patch File Download Diff File

45
mlprec/impl/mld_dmlprec_aply.f90
Unescape Escape View File

@ -784,6 +784,7 @@ contains
if (post) then if (post) then
if (.false.) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,& & dzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
@ -814,6 +815,26 @@ contains
& sweeps,work,info,init='Z') & sweeps,work,info,init='Z')
end if end if
else
!
! Apply the second smoother
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
end if
end if
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 POST smoother_apply') & a_err='Error during POST smoother_apply')
@ -871,6 +892,8 @@ contains
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post logical :: pre, post
real(psb_dpk_) :: tnrm
logical, parameter :: debug=.true.
character(len=20) :: name character(len=20) :: name
@ -933,7 +956,10 @@ contains
& a_err='Error during 2-PRE smoother_apply') & a_err='Error during 2-PRE smoother_apply')
goto 9999 goto 9999
end if end if
if (debug) then
tnrm = psb_genrm2(mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
write(0,*)' After first smoother y2l ',tnrm
end if
! !
! Compute the residual and call recursively ! Compute the residual and call recursively
@ -951,11 +977,20 @@ contains
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
if (debug) then
tnrm = psb_genrm2(mlprec_wrk(level)%vty,p%precv(level)%base_desc,info)
write(0,*)' Residual before restriction ',tnrm
end if
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level)%vty,& call psb_map_X2Y(done,mlprec_wrk(level)%vty,&
& dzero,mlprec_wrk(level + 1)%vx2l,& & dzero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work)
if (debug) then
tnrm = psb_genrm2(mlprec_wrk(level+1)%vx2l,p%precv(level+1)%base_desc,info)
write(0,*)' Output of restriction ',tnrm
end if
call mlprec_wrk(level + 1)%vy2l%zero() call mlprec_wrk(level + 1)%vy2l%zero()
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,&
@ -1070,6 +1105,7 @@ contains
real(psb_dpk_) :: l2_norm, delta, rtol=0.25 real(psb_dpk_) :: l2_norm, delta, rtol=0.25
real(psb_dpk_), allocatable :: temp_v(:) real(psb_dpk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
logical, parameter :: debug=.true.
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
@ -1105,8 +1141,8 @@ contains
goto 9999 goto 9999
end if end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info) delta = psb_genrm2(w, p%precv(level)%base_desc, info)
if (debug) write(0,*)'On entry delta ',delta
!Apply the preconditioner !Apply the preconditioner
call mlprec_wrk(level)%vy2l%set(dzero) call mlprec_wrk(level)%vy2l%set(dzero)
@ -1149,7 +1185,8 @@ contains
l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info) l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info)
iter = 0 iter = 0
if (debug) write(0,*) 'Alpha ',alpha,' l2_norm',l2_norm,&
& ' delta_old',delta_old,' tau',tau
if (l2_norm <= rtol*delta) then if (l2_norm <= rtol*delta) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), dzero, x, p%precv(level)%base_desc, info) call psb_geaxpby(alpha, d(idx), dzero, x, p%precv(level)%base_desc, info)

54
mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -115,25 +115,25 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
endif endif
if (sweeps == 0) then !!$ if (sweeps == 0) then
!!$
! !!$ !
! K^0 = I !!$ ! K^0 = I
! zero sweeps of any smoother is just the identity. !!$ ! zero sweeps of any smoother is just the identity.
! !!$ !
call psb_geaxpby(alpha,x,beta,y,desc_data,info) !!$ call psb_geaxpby(alpha,x,beta,y,desc_data,info)
!!$
else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then !!$ else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init !!$ ! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) !!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
if (info /= psb_success_) then !!$ if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,& !!$ call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in sub_aply Jacobi Sweeps = 1') !!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 !!$ goto 9999
endif !!$ endif
!!$
else if (sweeps >= 1) then !!$ else if (sweeps >= 1) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -196,14 +196,14 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
else !!$ else
!!$
info = psb_err_iarg_neg_ !!$ info = psb_err_iarg_neg_
call psb_errpush(info,name,& !!$ call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/)) !!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999 !!$ goto 9999
!!$
endif !!$ endif
if (n_col <= size(work)) then if (n_col <= size(work)) then

54
mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -115,25 +115,25 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
endif endif
if (sweeps == 0) then !!$ if (sweeps == 0) then
!!$
! !!$ !
! K^0 = I !!$ ! K^0 = I
! zero sweeps of any smoother is just the identity. !!$ ! zero sweeps of any smoother is just the identity.
! !!$ !
call psb_geaxpby(alpha,x,beta,y,desc_data,info) !!$ call psb_geaxpby(alpha,x,beta,y,desc_data,info)
!!$
else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then !!$ else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init !!$ ! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) !!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
if (info /= psb_success_) then !!$ if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,& !!$ call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in sub_aply Jacobi Sweeps = 1') !!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 !!$ goto 9999
endif !!$ endif
!!$
else if (sweeps >= 1) then !!$ else if (sweeps >= 1) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -196,14 +196,14 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
else !!$ else
!!$
info = psb_err_iarg_neg_ !!$ info = psb_err_iarg_neg_
call psb_errpush(info,name,& !!$ call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/)) !!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999 !!$ goto 9999
!!$
endif !!$ endif
if (n_col <= size(work)) then if (n_col <= size(work)) then

54
mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -115,25 +115,25 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
endif endif
if (sweeps == 0) then !!$ if (sweeps == 0) then
!!$
! !!$ !
! K^0 = I !!$ ! K^0 = I
! zero sweeps of any smoother is just the identity. !!$ ! zero sweeps of any smoother is just the identity.
! !!$ !
call psb_geaxpby(alpha,x,beta,y,desc_data,info) !!$ call psb_geaxpby(alpha,x,beta,y,desc_data,info)
!!$
else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then !!$ else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init !!$ ! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) !!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
if (info /= psb_success_) then !!$ if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,& !!$ call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in sub_aply Jacobi Sweeps = 1') !!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 !!$ goto 9999
endif !!$ endif
!!$
else if (sweeps >= 1) then !!$ else if (sweeps >= 1) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -196,14 +196,14 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
else !!$ else
!!$
info = psb_err_iarg_neg_ !!$ info = psb_err_iarg_neg_
call psb_errpush(info,name,& !!$ call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/)) !!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999 !!$ goto 9999
!!$
endif !!$ endif
if (n_col <= size(work)) then if (n_col <= size(work)) then

54
mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -115,25 +115,25 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if end if
endif endif
if (sweeps == 0) then !!$ if (sweeps == 0) then
!!$
! !!$ !
! K^0 = I !!$ ! K^0 = I
! zero sweeps of any smoother is just the identity. !!$ ! zero sweeps of any smoother is just the identity.
! !!$ !
call psb_geaxpby(alpha,x,beta,y,desc_data,info) !!$ call psb_geaxpby(alpha,x,beta,y,desc_data,info)
!!$
else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then !!$ else if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init !!$ ! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) !!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
!!$
if (info /= psb_success_) then !!$ if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,& !!$ call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in sub_aply Jacobi Sweeps = 1') !!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 !!$ goto 9999
endif !!$ endif
!!$
else if (sweeps >= 1) then !!$ else if (sweeps >= 1) then
! !
! !
! Apply multiple sweeps of a block-Jacobi solver ! Apply multiple sweeps of a block-Jacobi solver
@ -196,14 +196,14 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999 goto 9999
end if end if
else !!$ else
!!$
info = psb_err_iarg_neg_ !!$ info = psb_err_iarg_neg_
call psb_errpush(info,name,& !!$ call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/)) !!$ & i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999 !!$ goto 9999
!!$
endif !!$ endif
if (n_col <= size(work)) then if (n_col <= size(work)) then

Write Preview
Loading…
Cancel
Save