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

Use ASSOCIATE for wrk vectors. KCYCLE to be debugged.

Browse Source
stopcriterion
Salvatore Filippone 7 years ago
parent 08040c455b
commit 6f9a3c10d2
8 changed files with 2072 additions and 1932 deletions
Show Stats Download Patch File Download Diff File

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

@ -251,14 +251,18 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!
! At first iteration we must use the input BETA
!
beta_ = beta
level = 1
call psb_geaxpby(cone,x,czero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(cone,x,czero,vx2l,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
@ -276,17 +280,14 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = cone
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(cone,x,czero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-cone,p%precv(level)%base_a,y,&
& cone,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(cone,x,czero,vx2l,base_desc,info)
call psb_spmm(-cone,base_a,y,cone,vx2l,base_desc,info)
end do
!
@ -305,9 +306,9 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
end associate
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
@ -479,29 +480,33 @@ contains
goto 9999
end if
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
& vx2l,czero,vy2l,&
& base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vy2l,czero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& vy2l,czero,vty,&
& base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(cone,&
& p%wrk(level)%vtx,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,czero,vy2l,&
& base_desc, trans,&
& ione,work,info,init='Z')
end do
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,czero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
if (info /= psb_success_) then
@ -512,8 +517,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(cone,vx2l,&
& czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -531,8 +536,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -542,6 +547,7 @@ contains
end if
end associate
call psb_erractionrestore(err_act)
return
@ -597,6 +603,9 @@ contains
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level < nlev) then
!
! Apply the first smoother
@ -607,14 +616,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,czero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,czero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -628,20 +637,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(cone,vx2l,&
& czero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,cone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-cone,base_a,&
& vy2l,cone,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call psb_map_X2Y(cone,p%wrk(level)%vty,&
& czero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(cone,vty,&
& czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -650,8 +659,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(cone,vx2l,&
& czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -665,8 +674,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -676,14 +685,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,cone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(cone,p%wrk(level)%vty,&
& czero,p%wrk(level+1)%vx2l,&
call psb_geaxpby(cone,vx2l,&
& czero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,base_a,&
& vy2l,cone,vty,&
& base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(cone,vty,&
& czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -693,8 +702,8 @@ contains
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -707,12 +716,12 @@ contains
if (post) then
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,&
& cone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(cone,vx2l,&
& czero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,base_a,&
& vy2l,&
& cone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -726,14 +735,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,cone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,cone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -749,8 +758,8 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,czero,vy2l,&
& base_desc, trans,&
& sweeps,work,info)
else
@ -760,6 +769,7 @@ contains
& a_err='Invalid LEVEL vs NLEV')
goto 9999
end if
end associate
call psb_erractionrestore(err_act)
return
@ -824,14 +834,17 @@ contains
!K cycle
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,czero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else if (level < nlev) then
@ -839,14 +852,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,czero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,czero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -861,13 +874,13 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(cone,vx2l,&
& czero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,cone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-cone,base_a,&
& vy2l,cone,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -875,8 +888,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(cone,p%wrk(level)%vty,&
& czero,p%wrk(level + 1)%vx2l,&
call psb_map_X2Y(cone,vty,&
& czero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@ -910,8 +923,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -923,11 +936,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-cone,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& cone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(cone,vx2l,&
& czero,vty,&
& base_desc,info)
call psb_spmm(-cone,base_a,vy2l,&
& cone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -940,14 +953,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,cone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,cone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -964,7 +977,7 @@ contains
goto 9999
endif
end associate
call psb_erractionrestore(err_act)
return
@ -999,59 +1012,63 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(cone,p%wrk(level)%vx2l,czero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(cone,p%wrk(level)%vx2l,czero,w,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(cone,vx2l,czero,rhs,&
& base_desc,info)
call psb_geaxpby(cone,vx2l,czero,w,&
& base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
nc2l = base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
delta0 = psb_genrm2(w, base_desc, info)
!Apply the preconditioner
call p%wrk(level)%vy2l%zero()
call vy2l%zero()
idx=0
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(cone,p%wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(cone,vy2l,czero,d(idx),base_desc,info)
call psb_spmm(cone,p%precv(level)%base_a,d(idx),czero,v,p%precv(level)%base_desc,info)
call psb_spmm(cone,base_a,d(idx),czero,v,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1060,12 +1077,12 @@ contains
!FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(d(idx), w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info)
!GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(v, w, base_desc, info)
tau = psb_gedot(v, v, base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
@ -1074,26 +1091,26 @@ contains
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, cone, w, p%precv(level)%base_desc, info)
call psb_geaxpby(-alpha, v, cone, w, base_desc, info)
l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
l2_norm = psb_genrm2(w, base_desc, info)
iter = 0
if (l2_norm <= rtol*delta0) then
!Update solution x
call psb_geaxpby(alpha, d(idx), cone, x, p%precv(level)%base_desc, info)
call psb_geaxpby(alpha, d(idx), cone, x, base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(cone,w,czero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(cone,w,czero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(cone,p%wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(cone,vy2l,czero,d(idx),base_desc,info)
!Sparse matrix vector product
call psb_spmm(cone,p%precv(level)%base_a,d(idx),czero,v1,p%precv(level)%base_desc,info)
call psb_spmm(cone,base_a,d(idx),czero,v1,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1102,14 +1119,14 @@ contains
!tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau1= psb_gedot(d(idx), v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau1= psb_gedot(v1, v, base_desc, info)
tau2= psb_gedot(v1, v1, base_desc, info)
tau3= psb_gedot(v1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
@ -1119,20 +1136,20 @@ contains
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),cone,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx - 1),cone,x,base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),cone,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx),cone,x,base_desc,info)
endif
call psb_geaxpby(cone,x,czero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(cone,x,czero,vy2l,base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
call psb_gefree(v, 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(d(0), base_desc, info)
call psb_gefree(d(1), base_desc, info)
end associate
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then

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

@ -251,14 +251,18 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!
! At first iteration we must use the input BETA
!
beta_ = beta
level = 1
call psb_geaxpby(done,x,dzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(done,x,dzero,vx2l,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
@ -276,17 +280,14 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = done
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(done,x,dzero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-done,p%precv(level)%base_a,y,&
& done,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(done,x,dzero,vx2l,base_desc,info)
call psb_spmm(-done,base_a,y,done,vx2l,base_desc,info)
end do
!
@ -305,9 +306,9 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
end associate
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
@ -479,29 +480,33 @@ contains
goto 9999
end if
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
& vx2l,dzero,vy2l,&
& base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vy2l,dzero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& vy2l,dzero,vty,&
& base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(done,&
& p%wrk(level)%vtx,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,dzero,vy2l,&
& base_desc, trans,&
& ione,work,info,init='Z')
end do
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,dzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
if (info /= psb_success_) then
@ -512,8 +517,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(done,vx2l,&
& dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -531,8 +536,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -542,6 +547,7 @@ contains
end if
end associate
call psb_erractionrestore(err_act)
return
@ -597,6 +603,9 @@ contains
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level < nlev) then
!
! Apply the first smoother
@ -607,14 +616,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,dzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,dzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -628,20 +637,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(done,vx2l,&
& dzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,done,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call psb_map_X2Y(done,p%wrk(level)%vty,&
& dzero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -650,8 +659,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(done,vx2l,&
& dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -665,8 +674,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -676,14 +685,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,done,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(done,p%wrk(level)%vty,&
& dzero,p%wrk(level+1)%vx2l,&
call psb_geaxpby(done,vx2l,&
& dzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,&
& base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -693,8 +702,8 @@ contains
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -707,12 +716,12 @@ contains
if (post) then
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,&
& done,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(done,vx2l,&
& dzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,&
& done,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -726,14 +735,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,done,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,done,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -749,8 +758,8 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,dzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info)
else
@ -760,6 +769,7 @@ contains
& a_err='Invalid LEVEL vs NLEV')
goto 9999
end if
end associate
call psb_erractionrestore(err_act)
return
@ -824,14 +834,17 @@ contains
!K cycle
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,dzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else if (level < nlev) then
@ -839,14 +852,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,dzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,dzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -861,13 +874,13 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(done,vx2l,&
& dzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,done,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-done,base_a,&
& vy2l,done,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -875,8 +888,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(done,p%wrk(level)%vty,&
& dzero,p%wrk(level + 1)%vx2l,&
call psb_map_X2Y(done,vty,&
& dzero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@ -910,8 +923,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -923,11 +936,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-done,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& done,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(done,vx2l,&
& dzero,vty,&
& base_desc,info)
call psb_spmm(-done,base_a,vy2l,&
& done,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -940,14 +953,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,done,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,done,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -964,7 +977,7 @@ contains
goto 9999
endif
end associate
call psb_erractionrestore(err_act)
return
@ -999,59 +1012,63 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(done,p%wrk(level)%vx2l,dzero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(done,p%wrk(level)%vx2l,dzero,w,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(done,vx2l,dzero,rhs,&
& base_desc,info)
call psb_geaxpby(done,vx2l,dzero,w,&
& base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
nc2l = base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
delta0 = psb_genrm2(w, base_desc, info)
!Apply the preconditioner
call p%wrk(level)%vy2l%zero()
call vy2l%zero()
idx=0
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(done,p%wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(done,vy2l,dzero,d(idx),base_desc,info)
call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v,p%precv(level)%base_desc,info)
call psb_spmm(done,base_a,d(idx),dzero,v,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1060,12 +1077,12 @@ contains
!FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(d(idx), w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info)
!GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(v, w, base_desc, info)
tau = psb_gedot(v, v, base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
@ -1074,26 +1091,26 @@ contains
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, done, w, p%precv(level)%base_desc, info)
call psb_geaxpby(-alpha, v, done, w, base_desc, info)
l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
l2_norm = psb_genrm2(w, base_desc, info)
iter = 0
if (l2_norm <= rtol*delta0) then
!Update solution x
call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info)
call psb_geaxpby(alpha, d(idx), done, x, base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(done,w,dzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(done,w,dzero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(done,p%wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(done,vy2l,dzero,d(idx),base_desc,info)
!Sparse matrix vector product
call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v1,p%precv(level)%base_desc,info)
call psb_spmm(done,base_a,d(idx),dzero,v1,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1102,14 +1119,14 @@ contains
!tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau1= psb_gedot(d(idx), v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau1= psb_gedot(v1, v, base_desc, info)
tau2= psb_gedot(v1, v1, base_desc, info)
tau3= psb_gedot(v1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
@ -1119,20 +1136,20 @@ contains
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),done,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx - 1),done,x,base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx),done,x,base_desc,info)
endif
call psb_geaxpby(done,x,dzero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(done,x,dzero,vy2l,base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
call psb_gefree(v, 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(d(0), base_desc, info)
call psb_gefree(d(1), base_desc, info)
end associate
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then

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

@ -251,14 +251,18 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!
! At first iteration we must use the input BETA
!
beta_ = beta
level = 1
call psb_geaxpby(sone,x,szero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(sone,x,szero,vx2l,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
@ -276,17 +280,14 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = sone
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(sone,x,szero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-sone,p%precv(level)%base_a,y,&
& sone,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(sone,x,szero,vx2l,base_desc,info)
call psb_spmm(-sone,base_a,y,sone,vx2l,base_desc,info)
end do
!
@ -305,9 +306,9 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
end associate
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
@ -479,29 +480,33 @@ contains
goto 9999
end if
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
& vx2l,szero,vy2l,&
& base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vy2l,szero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& vy2l,szero,vty,&
& base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(sone,&
& p%wrk(level)%vtx,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,szero,vy2l,&
& base_desc, trans,&
& ione,work,info,init='Z')
end do
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,szero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
if (info /= psb_success_) then
@ -512,8 +517,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(sone,vx2l,&
& szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -531,8 +536,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -542,6 +547,7 @@ contains
end if
end associate
call psb_erractionrestore(err_act)
return
@ -597,6 +603,9 @@ contains
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level < nlev) then
!
! Apply the first smoother
@ -607,14 +616,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,szero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,szero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -628,20 +637,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(sone,vx2l,&
& szero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,sone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-sone,base_a,&
& vy2l,sone,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call psb_map_X2Y(sone,p%wrk(level)%vty,&
& szero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(sone,vty,&
& szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -650,8 +659,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(sone,vx2l,&
& szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -665,8 +674,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -676,14 +685,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,sone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(sone,p%wrk(level)%vty,&
& szero,p%wrk(level+1)%vx2l,&
call psb_geaxpby(sone,vx2l,&
& szero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,base_a,&
& vy2l,sone,vty,&
& base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(sone,vty,&
& szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -693,8 +702,8 @@ contains
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -707,12 +716,12 @@ contains
if (post) then
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,&
& sone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(sone,vx2l,&
& szero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,base_a,&
& vy2l,&
& sone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -726,14 +735,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,sone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,sone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -749,8 +758,8 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,szero,vy2l,&
& base_desc, trans,&
& sweeps,work,info)
else
@ -760,6 +769,7 @@ contains
& a_err='Invalid LEVEL vs NLEV')
goto 9999
end if
end associate
call psb_erractionrestore(err_act)
return
@ -824,14 +834,17 @@ contains
!K cycle
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,szero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else if (level < nlev) then
@ -839,14 +852,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,szero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,szero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -861,13 +874,13 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(sone,vx2l,&
& szero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,sone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-sone,base_a,&
& vy2l,sone,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -875,8 +888,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(sone,p%wrk(level)%vty,&
& szero,p%wrk(level + 1)%vx2l,&
call psb_map_X2Y(sone,vty,&
& szero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@ -910,8 +923,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -923,11 +936,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-sone,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& sone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(sone,vx2l,&
& szero,vty,&
& base_desc,info)
call psb_spmm(-sone,base_a,vy2l,&
& sone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -940,14 +953,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,sone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,sone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -964,7 +977,7 @@ contains
goto 9999
endif
end associate
call psb_erractionrestore(err_act)
return
@ -999,59 +1012,63 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(sone,p%wrk(level)%vx2l,szero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(sone,p%wrk(level)%vx2l,szero,w,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(sone,vx2l,szero,rhs,&
& base_desc,info)
call psb_geaxpby(sone,vx2l,szero,w,&
& base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
nc2l = base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
delta0 = psb_genrm2(w, base_desc, info)
!Apply the preconditioner
call p%wrk(level)%vy2l%zero()
call vy2l%zero()
idx=0
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(sone,p%wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(sone,vy2l,szero,d(idx),base_desc,info)
call psb_spmm(sone,p%precv(level)%base_a,d(idx),szero,v,p%precv(level)%base_desc,info)
call psb_spmm(sone,base_a,d(idx),szero,v,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1060,12 +1077,12 @@ contains
!FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(d(idx), w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info)
!GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(v, w, base_desc, info)
tau = psb_gedot(v, v, base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
@ -1074,26 +1091,26 @@ contains
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, sone, w, p%precv(level)%base_desc, info)
call psb_geaxpby(-alpha, v, sone, w, base_desc, info)
l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
l2_norm = psb_genrm2(w, base_desc, info)
iter = 0
if (l2_norm <= rtol*delta0) then
!Update solution x
call psb_geaxpby(alpha, d(idx), sone, x, p%precv(level)%base_desc, info)
call psb_geaxpby(alpha, d(idx), sone, x, base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(sone,w,szero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(sone,w,szero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(sone,p%wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(sone,vy2l,szero,d(idx),base_desc,info)
!Sparse matrix vector product
call psb_spmm(sone,p%precv(level)%base_a,d(idx),szero,v1,p%precv(level)%base_desc,info)
call psb_spmm(sone,base_a,d(idx),szero,v1,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1102,14 +1119,14 @@ contains
!tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau1= psb_gedot(d(idx), v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau1= psb_gedot(v1, v, base_desc, info)
tau2= psb_gedot(v1, v1, base_desc, info)
tau3= psb_gedot(v1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
@ -1119,20 +1136,20 @@ contains
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),sone,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx - 1),sone,x,base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),sone,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx),sone,x,base_desc,info)
endif
call psb_geaxpby(sone,x,szero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(sone,x,szero,vy2l,base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
call psb_gefree(v, 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(d(0), base_desc, info)
call psb_gefree(d(1), base_desc, info)
end associate
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then

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

@ -251,14 +251,18 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!
! At first iteration we must use the input BETA
!
beta_ = beta
level = 1
call psb_geaxpby(zone,x,zzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(zone,x,zzero,vx2l,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
@ -276,17 +280,14 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = zone
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(zone,x,zzero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-zone,p%precv(level)%base_a,y,&
& zone,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(zone,x,zzero,vx2l,base_desc,info)
call psb_spmm(-zone,base_a,y,zone,vx2l,base_desc,info)
end do
!
@ -305,9 +306,9 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,vy2l,beta_,y,base_desc,info)
end associate
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
@ -479,29 +480,33 @@ contains
goto 9999
end if
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
& vx2l,zzero,vy2l,&
& base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vy2l,zzero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& vy2l,zzero,vty,&
& base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(zone,&
& p%wrk(level)%vtx,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,zzero,vy2l,&
& base_desc, trans,&
& ione,work,info,init='Z')
end do
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,zzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
if (info /= psb_success_) then
@ -512,8 +517,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(zone,vx2l,&
& zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -531,8 +536,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -542,6 +547,7 @@ contains
end if
end associate
call psb_erractionrestore(err_act)
return
@ -597,6 +603,9 @@ contains
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level < nlev) then
!
! Apply the first smoother
@ -607,14 +616,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,zzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,zzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -628,20 +637,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(zone,vx2l,&
& zzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,zone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-zone,base_a,&
& vy2l,zone,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call psb_map_X2Y(zone,p%wrk(level)%vty,&
& zzero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(zone,vty,&
& zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -650,8 +659,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level+1)%vx2l,&
call psb_map_X2Y(zone,vx2l,&
& zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -665,8 +674,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -676,14 +685,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,zone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(zone,p%wrk(level)%vty,&
& zzero,p%wrk(level+1)%vx2l,&
call psb_geaxpby(zone,vx2l,&
& zzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,base_a,&
& vy2l,zone,vty,&
& base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(zone,vty,&
& zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -693,8 +702,8 @@ contains
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -707,12 +716,12 @@ contains
if (post) then
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,&
& zone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(zone,vx2l,&
& zzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,base_a,&
& vy2l,&
& zone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -726,14 +735,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,zone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,zone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -749,8 +758,8 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,zzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info)
else
@ -760,6 +769,7 @@ contains
& a_err='Invalid LEVEL vs NLEV')
goto 9999
end if
end associate
call psb_erractionrestore(err_act)
return
@ -824,14 +834,17 @@ contains
!K cycle
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
if (level == nlev) then
!
! Apply smoother
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,zzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else if (level < nlev) then
@ -839,14 +852,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,zzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vx2l,zzero,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -861,13 +874,13 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(zone,vx2l,&
& zzero,vty,&
& base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& p%wrk(level)%vy2l,zone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_spmm(-zone,base_a,&
& vy2l,zone,vty,&
& base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -875,8 +888,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(zone,p%wrk(level)%vty,&
& zzero,p%wrk(level + 1)%vx2l,&
call psb_map_X2Y(zone,vty,&
& zzero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@ -910,8 +923,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@ -923,11 +936,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-zone,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& zone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_geaxpby(zone,vx2l,&
& zzero,vty,&
& base_desc,info)
call psb_spmm(-zone,base_a,vy2l,&
& zone,vty,base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -940,14 +953,14 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,zone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& vty,zone,vy2l,&
& base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@ -964,7 +977,7 @@ contains
goto 9999
endif
end associate
call psb_erractionrestore(err_act)
return
@ -999,59 +1012,63 @@ contains
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
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,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc)
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
& base_desc,info,&
& scratch=.true.,mold=vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
& base_desc,info,&
& scratch=.true.,mold=vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(zone,p%wrk(level)%vx2l,zzero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(zone,p%wrk(level)%vx2l,zzero,w,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(zone,vx2l,zzero,rhs,&
& base_desc,info)
call psb_geaxpby(zone,vx2l,zzero,w,&
& base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
nc2l = base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
delta0 = psb_genrm2(w, base_desc, info)
!Apply the preconditioner
call p%wrk(level)%vy2l%zero()
call vy2l%zero()
idx=0
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(zone,p%wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(zone,vy2l,zzero,d(idx),base_desc,info)
call psb_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v,p%precv(level)%base_desc,info)
call psb_spmm(zone,base_a,d(idx),zzero,v,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1060,12 +1077,12 @@ contains
!FCG
if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(d(idx), w, base_desc, info)
tau = psb_gedot(d(idx), v, base_desc, info)
!GCR
else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
delta_old = psb_gedot(v, w, base_desc, info)
tau = psb_gedot(v, v, base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
@ -1074,26 +1091,26 @@ contains
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, zone, w, p%precv(level)%base_desc, info)
call psb_geaxpby(-alpha, v, zone, w, base_desc, info)
l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
l2_norm = psb_genrm2(w, base_desc, info)
iter = 0
if (l2_norm <= rtol*delta0) then
!Update solution x
call psb_geaxpby(alpha, d(idx), zone, x, p%precv(level)%base_desc, info)
call psb_geaxpby(alpha, d(idx), zone, x, base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(zone,w,zzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(zone,w,zzero,vx2l,base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(zone,p%wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(zone,vy2l,zzero,d(idx),base_desc,info)
!Sparse matrix vector product
call psb_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v1,p%precv(level)%base_desc,info)
call psb_spmm(zone,base_a,d(idx),zzero,v1,base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
@ -1102,14 +1119,14 @@ contains
!tau1, tau2, tau3, tau4
if (psb_toupper(trim(innersolv)) == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau1= psb_gedot(d(idx), v, base_desc, info)
tau2= psb_gedot(d(idx), v1, base_desc, info)
tau3= psb_gedot(d(idx), w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else if (psb_toupper(trim(innersolv)) == 'GCR') then
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau1= psb_gedot(v1, v, base_desc, info)
tau2= psb_gedot(v1, v1, base_desc, info)
tau3= psb_gedot(v1, w, base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
@ -1119,20 +1136,20 @@ contains
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),zone,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx - 1),zone,x,base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),zone,x,p%precv(level)%base_desc,info)
call psb_geaxpby(alpha,d(idx),zone,x,base_desc,info)
endif
call psb_geaxpby(zone,x,zzero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(zone,x,zzero,vy2l,base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
call psb_gefree(v, 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(d(0), base_desc, info)
call psb_gefree(d(1), base_desc, info)
end associate
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then

24
mlprec/mld_c_onelev_mod.f90
Unescape Escape View File

@ -552,11 +552,12 @@ contains
!
! Now for the ML application itself
!
! We have VTX/VTY/VX2L/VY2L
! VTX/VTY/VX2L/VY2L are stored explicitly
!
val = val + 4
!
! plus some additions for specific ML/cycles
! additions for specific ML/cycles
!
select case(lv%parms%ml_cycle)
case(mld_add_ml_,mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
@ -586,6 +587,19 @@ contains
info = psb_success_
nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vy2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vtx,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end subroutine c_base_onelev_allocate_wrk
@ -597,6 +611,10 @@ contains
!
integer(psb_ipk_) :: nwv
info = psb_success_
call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info)
end subroutine c_base_onelev_free_wrk
end module mld_c_onelev_mod

24
mlprec/mld_d_onelev_mod.f90
Unescape Escape View File

@ -552,11 +552,12 @@ contains
!
! Now for the ML application itself
!
! We have VTX/VTY/VX2L/VY2L
! VTX/VTY/VX2L/VY2L are stored explicitly
!
val = val + 4
!
! plus some additions for specific ML/cycles
! additions for specific ML/cycles
!
select case(lv%parms%ml_cycle)
case(mld_add_ml_,mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
@ -586,6 +587,19 @@ contains
info = psb_success_
nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vy2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vtx,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end subroutine d_base_onelev_allocate_wrk
@ -597,6 +611,10 @@ contains
!
integer(psb_ipk_) :: nwv
info = psb_success_
call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info)
end subroutine d_base_onelev_free_wrk
end module mld_d_onelev_mod

24
mlprec/mld_s_onelev_mod.f90
Unescape Escape View File

@ -552,11 +552,12 @@ contains
!
! Now for the ML application itself
!
! We have VTX/VTY/VX2L/VY2L
! VTX/VTY/VX2L/VY2L are stored explicitly
!
val = val + 4
!
! plus some additions for specific ML/cycles
! additions for specific ML/cycles
!
select case(lv%parms%ml_cycle)
case(mld_add_ml_,mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
@ -586,6 +587,19 @@ contains
info = psb_success_
nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vy2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vtx,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end subroutine s_base_onelev_allocate_wrk
@ -597,6 +611,10 @@ contains
!
integer(psb_ipk_) :: nwv
info = psb_success_
call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info)
end subroutine s_base_onelev_free_wrk
end module mld_s_onelev_mod

24
mlprec/mld_z_onelev_mod.f90
Unescape Escape View File

@ -552,11 +552,12 @@ contains
!
! Now for the ML application itself
!
! We have VTX/VTY/VX2L/VY2L
! VTX/VTY/VX2L/VY2L are stored explicitly
!
val = val + 4
!
! plus some additions for specific ML/cycles
! additions for specific ML/cycles
!
select case(lv%parms%ml_cycle)
case(mld_add_ml_,mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
@ -586,6 +587,19 @@ contains
info = psb_success_
nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vy2l,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vtx,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end subroutine z_base_onelev_allocate_wrk
@ -597,6 +611,10 @@ contains
!
integer(psb_ipk_) :: nwv
info = psb_success_
call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info)
end subroutine z_base_onelev_free_wrk
end module mld_z_onelev_mod

Write Preview
Loading…
Cancel
Save