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

version of K-cycle working, but to be investigated further.

Browse Source
stopcriterion
Salvatore Filippone 7 years ago
parent 2481fec23d
commit 68f5691a99
8 changed files with 196 additions and 148 deletions
Show Stats Download Patch File Download Diff File

81
mlprec/impl/mld_cmlprec_aply.f90
Unescape Escape View File

@ -847,7 +847,7 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -1017,7 +1017,7 @@ contains
!Other variables !Other variables
type(psb_c_vect_type) :: v, w, rhs, v1, x type(psb_c_vect_type) :: v, w, rhs, v1, x
type(psb_c_vect_type), dimension(0:1) :: d type(psb_c_vect_type) :: d0, d1
complex(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta complex(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_spk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm real(psb_spk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
@ -1025,37 +1025,50 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle' character(len=20) :: name = 'innerit_k_cycle'
if (size(p%precv(level)%wrk%wv)<7) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end if
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & v => p%precv(level)%wrk%wv(1), &
& w => p%precv(level)%wrk%wv(2),&
& rhs => p%precv(level)%wrk%wv(3))!, &
!!$ & v1 => p%precv(level)%wrk%wv(4), &
!!$ & x => p%precv(level)%wrk%wv(5), &
!!$ & d0 => p%precv(level)%wrk%wv(1), &
!!$ & d1 => p%precv(level)%wrk%wv(2))
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
call psb_geasb(rhs,& !!$ call psb_geasb(rhs,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(w,& !!$ call psb_geasb(w,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v,& !!$ call psb_geasb(v,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,& call psb_geasb(v1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
call psb_geasb(x,& call psb_geasb(x,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1) !Assemble d0 and d1
call psb_geasb(d(0),& call psb_geasb(d0,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),& call psb_geasb(d1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call x%zero() call x%zero()
! rhs=vx2l and w=rhs ! rhs=vx2l and w=rhs
@ -1080,9 +1093,9 @@ contains
idx=0 idx=0
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(cone,vy2l,czero,d(idx),base_desc,info) call psb_geaxpby(cone,vy2l,czero,d0,base_desc,info)
call psb_spmm(cone,base_a,d(idx),czero,v,base_desc,info) call psb_spmm(cone,base_a,d0,czero,v,base_desc,info)
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')
@ -1091,8 +1104,8 @@ contains
!FCG !FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, base_desc, info) delta_old = psb_gedot(d0, w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info) tau = psb_gedot(d0, v, base_desc, info)
!GCR !GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, base_desc, info) delta_old = psb_gedot(v, w, base_desc, info)
@ -1112,7 +1125,7 @@ contains
if (l2_norm <= rtol*delta0) then if (l2_norm <= rtol*delta0) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), cone, x, base_desc, info) call psb_geaxpby(alpha, d0, cone, x, base_desc, info)
else else
iter = iter + 1 iter = iter + 1
idx=mod(iter,2) idx=mod(iter,2)
@ -1120,11 +1133,11 @@ contains
!Apply preconditioner !Apply preconditioner
call psb_geaxpby(cone,w,czero,vx2l,base_desc,info) call psb_geaxpby(cone,w,czero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(cone,vy2l,czero,d(idx),base_desc,info) call psb_geaxpby(cone,vy2l,czero,d1,base_desc,info)
!Sparse matrix vector product !Sparse matrix vector product
call psb_spmm(cone,base_a,d(idx),czero,v1,base_desc,info) call psb_spmm(cone,base_a,d1,czero,v1,base_desc,info)
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')
@ -1133,9 +1146,9 @@ contains
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, base_desc, info) tau1= psb_gedot(d1, v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info) tau2= psb_gedot(d1, v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info) tau3= psb_gedot(d1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, base_desc, info) tau1= psb_gedot(v1, v, base_desc, info)
@ -1150,19 +1163,20 @@ contains
!Update solution !Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4) alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),cone,x,base_desc,info) call psb_geaxpby(alpha,d0,cone,x,base_desc,info)
alpha=tau3/tau4 alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),cone,x,base_desc,info) call psb_geaxpby(alpha,d1,cone,x,base_desc,info)
endif endif
call psb_geaxpby(cone,x,czero,vy2l,base_desc,info) call psb_geaxpby(cone,x,czero,vy2l,base_desc,info)
!Free vectors !Free vectors
call psb_gefree(v, base_desc, info) !!$ call psb_gefree(v, base_desc, info)
!!$ call psb_gefree(w, base_desc, info)
!!$ call psb_gefree(rhs, base_desc, info)
call psb_gefree(v1, base_desc, info) call psb_gefree(v1, base_desc, info)
call psb_gefree(w, base_desc, info)
call psb_gefree(x, base_desc, info) call psb_gefree(x, base_desc, info)
call psb_gefree(d(0), base_desc, info) call psb_gefree(d0, base_desc, info)
call psb_gefree(d(1), base_desc, info) call psb_gefree(d1, base_desc, info)
end associate end associate
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1176,9 +1190,6 @@ contains
end subroutine mld_cmlprec_aply_vect end subroutine mld_cmlprec_aply_vect
! !
! Old routine for arrays instead of psb_X_vector. To be deleted eventually. ! Old routine for arrays instead of psb_X_vector. To be deleted eventually.
! !

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

@ -847,7 +847,7 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -1017,7 +1017,7 @@ contains
!Other variables !Other variables
type(psb_d_vect_type) :: v, w, rhs, v1, x type(psb_d_vect_type) :: v, w, rhs, v1, x
type(psb_d_vect_type), dimension(0:1) :: d type(psb_d_vect_type) :: d0, d1
real(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta real(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
@ -1025,37 +1025,50 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle' character(len=20) :: name = 'innerit_k_cycle'
if (size(p%precv(level)%wrk%wv)<7) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end if
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & v => p%precv(level)%wrk%wv(1), &
& w => p%precv(level)%wrk%wv(2),&
& rhs => p%precv(level)%wrk%wv(3))!, &
!!$ & v1 => p%precv(level)%wrk%wv(4), &
!!$ & x => p%precv(level)%wrk%wv(5), &
!!$ & d0 => p%precv(level)%wrk%wv(1), &
!!$ & d1 => p%precv(level)%wrk%wv(2))
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
call psb_geasb(rhs,& !!$ call psb_geasb(rhs,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(w,& !!$ call psb_geasb(w,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v,& !!$ call psb_geasb(v,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,& call psb_geasb(v1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
call psb_geasb(x,& call psb_geasb(x,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1) !Assemble d0 and d1
call psb_geasb(d(0),& call psb_geasb(d0,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),& call psb_geasb(d1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call x%zero() call x%zero()
! rhs=vx2l and w=rhs ! rhs=vx2l and w=rhs
@ -1080,9 +1093,9 @@ contains
idx=0 idx=0
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(done,vy2l,dzero,d(idx),base_desc,info) call psb_geaxpby(done,vy2l,dzero,d0,base_desc,info)
call psb_spmm(done,base_a,d(idx),dzero,v,base_desc,info) call psb_spmm(done,base_a,d0,dzero,v,base_desc,info)
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')
@ -1091,8 +1104,8 @@ contains
!FCG !FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, base_desc, info) delta_old = psb_gedot(d0, w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info) tau = psb_gedot(d0, v, base_desc, info)
!GCR !GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, base_desc, info) delta_old = psb_gedot(v, w, base_desc, info)
@ -1112,7 +1125,7 @@ contains
if (l2_norm <= rtol*delta0) then if (l2_norm <= rtol*delta0) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), done, x, base_desc, info) call psb_geaxpby(alpha, d0, done, x, base_desc, info)
else else
iter = iter + 1 iter = iter + 1
idx=mod(iter,2) idx=mod(iter,2)
@ -1120,11 +1133,11 @@ contains
!Apply preconditioner !Apply preconditioner
call psb_geaxpby(done,w,dzero,vx2l,base_desc,info) call psb_geaxpby(done,w,dzero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(done,vy2l,dzero,d(idx),base_desc,info) call psb_geaxpby(done,vy2l,dzero,d1,base_desc,info)
!Sparse matrix vector product !Sparse matrix vector product
call psb_spmm(done,base_a,d(idx),dzero,v1,base_desc,info) call psb_spmm(done,base_a,d1,dzero,v1,base_desc,info)
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')
@ -1133,9 +1146,9 @@ contains
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, base_desc, info) tau1= psb_gedot(d1, v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info) tau2= psb_gedot(d1, v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info) tau3= psb_gedot(d1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, base_desc, info) tau1= psb_gedot(v1, v, base_desc, info)
@ -1150,19 +1163,20 @@ contains
!Update solution !Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4) alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),done,x,base_desc,info) call psb_geaxpby(alpha,d0,done,x,base_desc,info)
alpha=tau3/tau4 alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),done,x,base_desc,info) call psb_geaxpby(alpha,d1,done,x,base_desc,info)
endif endif
call psb_geaxpby(done,x,dzero,vy2l,base_desc,info) call psb_geaxpby(done,x,dzero,vy2l,base_desc,info)
!Free vectors !Free vectors
call psb_gefree(v, base_desc, info) !!$ call psb_gefree(v, base_desc, info)
!!$ call psb_gefree(w, base_desc, info)
!!$ call psb_gefree(rhs, base_desc, info)
call psb_gefree(v1, base_desc, info) call psb_gefree(v1, base_desc, info)
call psb_gefree(w, base_desc, info)
call psb_gefree(x, base_desc, info) call psb_gefree(x, base_desc, info)
call psb_gefree(d(0), base_desc, info) call psb_gefree(d0, base_desc, info)
call psb_gefree(d(1), base_desc, info) call psb_gefree(d1, base_desc, info)
end associate end associate
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1176,9 +1190,6 @@ contains
end subroutine mld_dmlprec_aply_vect end subroutine mld_dmlprec_aply_vect
! !
! Old routine for arrays instead of psb_X_vector. To be deleted eventually. ! Old routine for arrays instead of psb_X_vector. To be deleted eventually.
! !

81
mlprec/impl/mld_smlprec_aply.f90
Unescape Escape View File

@ -847,7 +847,7 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -1017,7 +1017,7 @@ contains
!Other variables !Other variables
type(psb_s_vect_type) :: v, w, rhs, v1, x type(psb_s_vect_type) :: v, w, rhs, v1, x
type(psb_s_vect_type), dimension(0:1) :: d type(psb_s_vect_type) :: d0, d1
real(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta real(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_spk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm real(psb_spk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
@ -1025,37 +1025,50 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle' character(len=20) :: name = 'innerit_k_cycle'
if (size(p%precv(level)%wrk%wv)<7) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end if
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & v => p%precv(level)%wrk%wv(1), &
& w => p%precv(level)%wrk%wv(2),&
& rhs => p%precv(level)%wrk%wv(3))!, &
!!$ & v1 => p%precv(level)%wrk%wv(4), &
!!$ & x => p%precv(level)%wrk%wv(5), &
!!$ & d0 => p%precv(level)%wrk%wv(1), &
!!$ & d1 => p%precv(level)%wrk%wv(2))
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
call psb_geasb(rhs,& !!$ call psb_geasb(rhs,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(w,& !!$ call psb_geasb(w,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v,& !!$ call psb_geasb(v,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,& call psb_geasb(v1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
call psb_geasb(x,& call psb_geasb(x,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1) !Assemble d0 and d1
call psb_geasb(d(0),& call psb_geasb(d0,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),& call psb_geasb(d1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call x%zero() call x%zero()
! rhs=vx2l and w=rhs ! rhs=vx2l and w=rhs
@ -1080,9 +1093,9 @@ contains
idx=0 idx=0
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(sone,vy2l,szero,d(idx),base_desc,info) call psb_geaxpby(sone,vy2l,szero,d0,base_desc,info)
call psb_spmm(sone,base_a,d(idx),szero,v,base_desc,info) call psb_spmm(sone,base_a,d0,szero,v,base_desc,info)
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')
@ -1091,8 +1104,8 @@ contains
!FCG !FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, base_desc, info) delta_old = psb_gedot(d0, w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info) tau = psb_gedot(d0, v, base_desc, info)
!GCR !GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, base_desc, info) delta_old = psb_gedot(v, w, base_desc, info)
@ -1112,7 +1125,7 @@ contains
if (l2_norm <= rtol*delta0) then if (l2_norm <= rtol*delta0) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), sone, x, base_desc, info) call psb_geaxpby(alpha, d0, sone, x, base_desc, info)
else else
iter = iter + 1 iter = iter + 1
idx=mod(iter,2) idx=mod(iter,2)
@ -1120,11 +1133,11 @@ contains
!Apply preconditioner !Apply preconditioner
call psb_geaxpby(sone,w,szero,vx2l,base_desc,info) call psb_geaxpby(sone,w,szero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(sone,vy2l,szero,d(idx),base_desc,info) call psb_geaxpby(sone,vy2l,szero,d1,base_desc,info)
!Sparse matrix vector product !Sparse matrix vector product
call psb_spmm(sone,base_a,d(idx),szero,v1,base_desc,info) call psb_spmm(sone,base_a,d1,szero,v1,base_desc,info)
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')
@ -1133,9 +1146,9 @@ contains
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, base_desc, info) tau1= psb_gedot(d1, v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info) tau2= psb_gedot(d1, v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info) tau3= psb_gedot(d1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, base_desc, info) tau1= psb_gedot(v1, v, base_desc, info)
@ -1150,19 +1163,20 @@ contains
!Update solution !Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4) alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),sone,x,base_desc,info) call psb_geaxpby(alpha,d0,sone,x,base_desc,info)
alpha=tau3/tau4 alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),sone,x,base_desc,info) call psb_geaxpby(alpha,d1,sone,x,base_desc,info)
endif endif
call psb_geaxpby(sone,x,szero,vy2l,base_desc,info) call psb_geaxpby(sone,x,szero,vy2l,base_desc,info)
!Free vectors !Free vectors
call psb_gefree(v, base_desc, info) !!$ call psb_gefree(v, base_desc, info)
!!$ call psb_gefree(w, base_desc, info)
!!$ call psb_gefree(rhs, base_desc, info)
call psb_gefree(v1, base_desc, info) call psb_gefree(v1, base_desc, info)
call psb_gefree(w, base_desc, info)
call psb_gefree(x, base_desc, info) call psb_gefree(x, base_desc, info)
call psb_gefree(d(0), base_desc, info) call psb_gefree(d0, base_desc, info)
call psb_gefree(d(1), base_desc, info) call psb_gefree(d1, base_desc, info)
end associate end associate
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1176,9 +1190,6 @@ contains
end subroutine mld_smlprec_aply_vect end subroutine mld_smlprec_aply_vect
! !
! Old routine for arrays instead of psb_X_vector. To be deleted eventually. ! Old routine for arrays instead of psb_X_vector. To be deleted eventually.
! !

81
mlprec/impl/mld_zmlprec_aply.f90
Unescape Escape View File

@ -847,7 +847,7 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & wv => p%precv(level)%wrk%wv(8:))
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -1017,7 +1017,7 @@ contains
!Other variables !Other variables
type(psb_z_vect_type) :: v, w, rhs, v1, x type(psb_z_vect_type) :: v, w, rhs, v1, x
type(psb_z_vect_type), dimension(0:1) :: d type(psb_z_vect_type) :: d0, d1
complex(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta complex(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
@ -1025,37 +1025,50 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle' character(len=20) :: name = 'innerit_k_cycle'
if (size(p%precv(level)%wrk%wv)<7) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='invalid wv size')
goto 9999
end if
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv) & v => p%precv(level)%wrk%wv(1), &
& w => p%precv(level)%wrk%wv(2),&
& rhs => p%precv(level)%wrk%wv(3))!, &
!!$ & v1 => p%precv(level)%wrk%wv(4), &
!!$ & x => p%precv(level)%wrk%wv(5), &
!!$ & d0 => p%precv(level)%wrk%wv(1), &
!!$ & d1 => p%precv(level)%wrk%wv(2))
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
call psb_geasb(rhs,& !!$ call psb_geasb(rhs,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(w,& !!$ call psb_geasb(w,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v,& !!$ call psb_geasb(v,&
& base_desc,info,& !!$ & base_desc,info,&
& scratch=.true.,mold=vx2l%v) !!$ & scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,& call psb_geasb(v1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
call psb_geasb(x,& call psb_geasb(x,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vx2l%v) & scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1) !Assemble d0 and d1
call psb_geasb(d(0),& call psb_geasb(d0,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),& call psb_geasb(d1,&
& base_desc,info,& & base_desc,info,&
& scratch=.true.,mold=vy2l%v) & scratch=.true.,mold=vy2l%v)
call x%zero() call x%zero()
! rhs=vx2l and w=rhs ! rhs=vx2l and w=rhs
@ -1080,9 +1093,9 @@ contains
idx=0 idx=0
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(zone,vy2l,zzero,d(idx),base_desc,info) call psb_geaxpby(zone,vy2l,zzero,d0,base_desc,info)
call psb_spmm(zone,base_a,d(idx),zzero,v,base_desc,info) call psb_spmm(zone,base_a,d0,zzero,v,base_desc,info)
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')
@ -1091,8 +1104,8 @@ contains
!FCG !FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, base_desc, info) delta_old = psb_gedot(d0, w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info) tau = psb_gedot(d0, v, base_desc, info)
!GCR !GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, base_desc, info) delta_old = psb_gedot(v, w, base_desc, info)
@ -1112,7 +1125,7 @@ contains
if (l2_norm <= rtol*delta0) then if (l2_norm <= rtol*delta0) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), zone, x, base_desc, info) call psb_geaxpby(alpha, d0, zone, x, base_desc, info)
else else
iter = iter + 1 iter = iter + 1
idx=mod(iter,2) idx=mod(iter,2)
@ -1120,11 +1133,11 @@ contains
!Apply preconditioner !Apply preconditioner
call psb_geaxpby(zone,w,zzero,vx2l,base_desc,info) call psb_geaxpby(zone,w,zzero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info) call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(zone,vy2l,zzero,d(idx),base_desc,info) call psb_geaxpby(zone,vy2l,zzero,d1,base_desc,info)
!Sparse matrix vector product !Sparse matrix vector product
call psb_spmm(zone,base_a,d(idx),zzero,v1,base_desc,info) call psb_spmm(zone,base_a,d1,zzero,v1,base_desc,info)
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')
@ -1133,9 +1146,9 @@ contains
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, base_desc, info) tau1= psb_gedot(d1, v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info) tau2= psb_gedot(d1, v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info) tau3= psb_gedot(d1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, base_desc, info) tau1= psb_gedot(v1, v, base_desc, info)
@ -1150,19 +1163,20 @@ contains
!Update solution !Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4) alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),zone,x,base_desc,info) call psb_geaxpby(alpha,d0,zone,x,base_desc,info)
alpha=tau3/tau4 alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),zone,x,base_desc,info) call psb_geaxpby(alpha,d1,zone,x,base_desc,info)
endif endif
call psb_geaxpby(zone,x,zzero,vy2l,base_desc,info) call psb_geaxpby(zone,x,zzero,vy2l,base_desc,info)
!Free vectors !Free vectors
call psb_gefree(v, base_desc, info) !!$ call psb_gefree(v, base_desc, info)
!!$ call psb_gefree(w, base_desc, info)
!!$ call psb_gefree(rhs, base_desc, info)
call psb_gefree(v1, base_desc, info) call psb_gefree(v1, base_desc, info)
call psb_gefree(w, base_desc, info)
call psb_gefree(x, base_desc, info) call psb_gefree(x, base_desc, info)
call psb_gefree(d(0), base_desc, info) call psb_gefree(d0, base_desc, info)
call psb_gefree(d(1), base_desc, info) call psb_gefree(d1, base_desc, info)
end associate end associate
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -1176,9 +1190,6 @@ contains
end subroutine mld_zmlprec_aply_vect end subroutine mld_zmlprec_aply_vect
! !
! Old routine for arrays instead of psb_X_vector. To be deleted eventually. ! Old routine for arrays instead of psb_X_vector. To be deleted eventually.
! !

5
mlprec/mld_c_onelev_mod.f90
Unescape Escape View File

@ -573,9 +573,10 @@ contains
case(mld_kcycle_ml_, mld_kcyclesym_ml_) case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !
! We need 7 in inneritkcycle, but we can reuse vtx ! We need 7 in inneritkcycle.
! Can we reuse vtx?
! !
val = val + 6 val = val + 7
case default case default
! Need a better error signaling ? ! Need a better error signaling ?

5
mlprec/mld_d_onelev_mod.f90
Unescape Escape View File

@ -573,9 +573,10 @@ contains
case(mld_kcycle_ml_, mld_kcyclesym_ml_) case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !
! We need 7 in inneritkcycle, but we can reuse vtx ! We need 7 in inneritkcycle.
! Can we reuse vtx?
! !
val = val + 6 val = val + 7
case default case default
! Need a better error signaling ? ! Need a better error signaling ?

5
mlprec/mld_s_onelev_mod.f90
Unescape Escape View File

@ -573,9 +573,10 @@ contains
case(mld_kcycle_ml_, mld_kcyclesym_ml_) case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !
! We need 7 in inneritkcycle, but we can reuse vtx ! We need 7 in inneritkcycle.
! Can we reuse vtx?
! !
val = val + 6 val = val + 7
case default case default
! Need a better error signaling ? ! Need a better error signaling ?

5
mlprec/mld_z_onelev_mod.f90
Unescape Escape View File

@ -573,9 +573,10 @@ contains
case(mld_kcycle_ml_, mld_kcyclesym_ml_) case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !
! We need 7 in inneritkcycle, but we can reuse vtx ! We need 7 in inneritkcycle.
! Can we reuse vtx?
! !
val = val + 6 val = val + 7
case default case default
! Need a better error signaling ? ! Need a better error signaling ?

Write Preview
Loading…
Cancel
Save