From 6f9a3c10d246b66ac426f081586eda073f0dc271 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 8 Dec 2017 15:17:31 +0000 Subject: [PATCH] Use ASSOCIATE for wrk vectors. KCYCLE to be debugged. --- mlprec/impl/mld_cmlprec_aply.f90 | 977 ++++++++++++++++--------------- mlprec/impl/mld_dmlprec_aply.f90 | 977 ++++++++++++++++--------------- mlprec/impl/mld_smlprec_aply.f90 | 977 ++++++++++++++++--------------- mlprec/impl/mld_zmlprec_aply.f90 | 977 ++++++++++++++++--------------- mlprec/mld_c_onelev_mod.f90 | 24 +- mlprec/mld_d_onelev_mod.f90 | 24 +- mlprec/mld_s_onelev_mod.f90 | 24 +- mlprec/mld_z_onelev_mod.f90 | 24 +- 8 files changed, 2072 insertions(+), 1932 deletions(-) diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index 4378f3d6..11629f69 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -251,20 +251,50 @@ 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 - ! - ! 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy') - goto 9999 - end if - - do isweep = 1, p%outer_sweeps - 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 + + + 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 + end if + + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call p%wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + 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,vx2l,base_desc,info) + call psb_spmm(-cone,base_a,y,cone,vx2l,base_desc,info) + end do + + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + ! ! With the current implementation, y2l is zeroed internally at first smoother. ! call p%wrk(level)%vy2l%zero() @@ -276,38 +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) - ! 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) - end do - - ! - ! If outer_sweeps == 1 we have just skipped the loop, and it's - ! equivalent to a single application. - ! - - ! - ! With the current implementation, y2l is zeroed internally at first smoother. - ! call p%wrk(level)%vy2l%zero() - ! - call inner_ml_aply(level,p,trans_,work,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,70 +480,75 @@ contains goto 9999 end if - 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) + 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,& + & 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,& + & vy2l,czero,vty,& + & base_desc, trans,& + & ione,work,info,init='Z') + + call p%precv(level)%sm2a%apply(cone,& + & vty,czero,vy2l,& + & base_desc, trans,& + & ione,work,info,init='Z') + end do - sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) - do k=1, sweeps + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(cone,& - & p%wrk(level)%vy2l,czero,p%wrk(level)%vtx,& - & p%precv(level)%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,& - & 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - ! Apply the restriction - call psb_map_X2Y(cone,p%wrk(level)%vx2l,& - & czero,p%wrk(level+1)%vx2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 + & vx2l,czero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') end if - - call inner_ml_aply(level+1,p,trans,work,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') + & a_err='Error during ADD smoother_apply') goto 9999 end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,& - & cone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + if (level < nlev) then + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + ! + ! Apply the prolongator + ! + 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,& + & a_err='Error during prolongation') + goto 9999 + end if - end if + end if + end associate + call psb_erractionrestore(err_act) return @@ -597,170 +603,174 @@ 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')) - if (level < nlev) then - ! - ! Apply the first smoother - ! The residual has been prepared before the recursive call. - ! + 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 + ! The residual has been prepared before the recursive call. + ! - if (pre) then - 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,& - & sweeps,work,info,init='Z') + if (pre) then + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(cone,& + & 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,& + & vx2l,czero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during PRE smoother_apply') + goto 9999 + end if + endif + ! + ! Compute the residual and call recursively + ! + if (pre) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + 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,& + & a_err='Error during restriction') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - endif - ! - ! 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) - - 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_) 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if - else - ! Shortcut: just transfer x2l. - call psb_map_X2Y(cone,p%wrk(level)%vx2l,& - & czero,p%wrk(level+1)%vx2l,& + ! Shortcut: just transfer x2l. + 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,& + & a_err='Error during restriction') + goto 9999 + end if + endif + + call inner_ml_aply(level+1,p,trans,work,info) + + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during restriction') + & a_err='Error during prolongation') goto 9999 end if - endif - call inner_ml_aply(level+1,p,trans,work,info) + if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then + + 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,& + & a_err='Error during W-cycle restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + + 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 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during W recusion/prolongation') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,& - & cone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + endif - 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W-cycle restriction') - goto 9999 - end if - - 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,& - & p%precv(level+1)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W recusion/prolongation') - goto 9999 - end if - - endif - - - 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,& - & 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 - ! - ! Apply the second smoother - ! - if (trans == 'N') then - sweeps = p%precv(level)%parms%sweeps_post - if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& - & p%wrk(level)%vty,cone,p%wrk(level)%vy2l,& - & p%precv(level)%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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if - - endif - - else if (level == nlev) 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,& - & sweeps,work,info) + if (post) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + + ! + ! Apply the second smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& + & 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,& + & vty,cone,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if - else + endif - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - end if + else if (level == nlev) then + + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(cone,& + & vx2l,czero,vy2l,& + & base_desc, trans,& + & sweeps,work,info) + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') + goto 9999 + end if + end associate + call psb_erractionrestore(err_act) return @@ -824,147 +834,150 @@ contains !K cycle - 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,& - & sweeps,work,info,init='Z') - - else if (level < nlev) then - - if (trans == 'N') then + 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,& - & 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 - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during 2-PRE smoother_apply') - goto 9999 - end if - + else if (level < nlev) then - ! - ! Compute the residual and call recursively - ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(cone,& + & 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,& + & vx2l,czero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if - 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_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during 2-PRE smoother_apply') + goto 9999 + end if - ! Apply the restriction - call psb_map_X2Y(cone,p%wrk(level)%vty,& - & czero,p%wrk(level + 1)%vx2l,& - & p%precv(level + 1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if + ! + ! Compute the residual and call recursively + ! - !Set the preconditioner + call psb_geaxpby(cone,vx2l,& + & czero,vty,& + & base_desc,info) - if (level <= nlev - 2 ) then - if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then - call mld_cinneritkcycle(p, level + 1, trans, work, 'FCG') - elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then - call mld_cinneritkcycle(p, level + 1, trans, work, 'GCR') - else + 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='Bad value for ml_cycle') + & a_err='Error during residue') goto 9999 + end if + + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + !Set the preconditioner + + if (level <= nlev - 2 ) then + if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then + call mld_cinneritkcycle(p, level + 1, trans, work, 'FCG') + elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then + call mld_cinneritkcycle(p, level + 1, trans, work, 'GCR') + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Bad value for ml_cycle') + goto 9999 + endif + else + call inner_ml_aply(level + 1 ,p,trans,work,info) endif - else - call inner_ml_aply(level + 1 ,p,trans,work,info) - endif - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,& - & cone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during prolongation') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - ! - ! 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,& - & 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 - ! - ! Apply the smoother - ! - 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,& - & sweeps,work,info,init='Z') + ! + ! Compute the residual + ! + 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,& + & a_err='Error during residue') + goto 9999 + end if + ! + ! Apply the smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& + & 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,& + & vty,cone,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') goto 9999 - end if - else - - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - - endif + endif + end associate call psb_erractionrestore(err_act) return @@ -998,141 +1011,145 @@ contains complex(psb_spk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx character(len=20) :: name = 'innerit_k_cycle' - - !Assemble rhs, w, v, v1, x - - call psb_geasb(rhs,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(w,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v1,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(x,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - - if (psb_errstatus_fatal()) then - nc2l = p%precv(level)%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) - - !Apply the preconditioner - call p%wrk(level)%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_spmm(cone,p%precv(level)%base_a,d(idx),czero,v,p%precv(level)%base_desc,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if - - !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) - !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) - else - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid inner solver') - goto 9999 - endif - alpha = delta_old/tau - !Update residual w - call psb_geaxpby(-alpha, v, cone, w, p%precv(level)%base_desc, info) + 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,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(w,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v1,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(x,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + call psb_geasb(d(1),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + + + call x%zero() + + ! rhs=vx2l and w=rhs + 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 = 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 - l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) - iter = 0 + delta0 = psb_genrm2(w, base_desc, info) - if (l2_norm <= rtol*delta0) then - !Update solution x - call psb_geaxpby(alpha, d(idx), cone, x, p%precv(level)%base_desc, info) - else - iter = iter + 1 - idx=mod(iter,2) + !Apply the preconditioner + call vy2l%zero() - !Apply preconditioner - call psb_geaxpby(cone,w,czero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info) + 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) - !Sparse matrix vector product + call psb_geaxpby(cone,vy2l,czero,d(idx),base_desc,info) - 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,v,base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during residue') goto 9999 end if - !tau1, tau2, tau3, tau4 + !FCG 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) - tau4= tau2 - (tau1*tau1)/tau + 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 - 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) - tau4= tau2 - (tau1*tau1)/tau + 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') goto 9999 endif - !Update solution - alpha=alpha-(tau1*tau3)/(tau*tau4) - call psb_geaxpby(alpha,d(idx - 1),cone,x,p%precv(level)%base_desc,info) - alpha=tau3/tau4 - call psb_geaxpby(alpha,d(idx),cone,x,p%precv(level)%base_desc,info) - endif + alpha = delta_old/tau + !Update residual w + call psb_geaxpby(-alpha, v, cone, w, 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, base_desc, info) + else + iter = iter + 1 + idx=mod(iter,2) - call psb_geaxpby(cone,x,czero,p%wrk(level)%vy2l,p%precv(level)%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) + !Apply preconditioner + call psb_geaxpby(cone,w,czero,vx2l,base_desc,info) + call inner_ml_aply(level,p,trans,work,info) + call psb_geaxpby(cone,vy2l,czero,d(idx),base_desc,info) + + !Sparse matrix vector product + + 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') + goto 9999 + end if + + !tau1, tau2, tau3, tau4 + if (psb_toupper(trim(innersolv)) == 'FCG') then + 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, 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,& + & a_err='Invalid inner solver') + goto 9999 + endif + + !Update solution + alpha=alpha-(tau1*tau3)/(tau*tau4) + call psb_geaxpby(alpha,d(idx - 1),cone,x,base_desc,info) + alpha=tau3/tau4 + call psb_geaxpby(alpha,d(idx),cone,x,base_desc,info) + endif + call psb_geaxpby(cone,x,czero,vy2l,base_desc,info) + !Free vectors + 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 diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index 3a0010cd..4a046d80 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -251,20 +251,50 @@ 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 - ! - ! 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy') - goto 9999 - end if - - do isweep = 1, p%outer_sweeps - 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 + + + 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 + end if + + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call p%wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + 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,vx2l,base_desc,info) + call psb_spmm(-done,base_a,y,done,vx2l,base_desc,info) + end do + + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + ! ! With the current implementation, y2l is zeroed internally at first smoother. ! call p%wrk(level)%vy2l%zero() @@ -276,38 +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) - ! 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) - end do - - ! - ! If outer_sweeps == 1 we have just skipped the loop, and it's - ! equivalent to a single application. - ! - - ! - ! With the current implementation, y2l is zeroed internally at first smoother. - ! call p%wrk(level)%vy2l%zero() - ! - call inner_ml_aply(level,p,trans_,work,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,70 +480,75 @@ contains goto 9999 end if - 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) + 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,& + & 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,& + & vy2l,dzero,vty,& + & base_desc, trans,& + & ione,work,info,init='Z') + + call p%precv(level)%sm2a%apply(done,& + & vty,dzero,vy2l,& + & base_desc, trans,& + & ione,work,info,init='Z') + end do - sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) - do k=1, sweeps + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(done,& - & p%wrk(level)%vy2l,dzero,p%wrk(level)%vtx,& - & p%precv(level)%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,& - & 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - ! Apply the restriction - call psb_map_X2Y(done,p%wrk(level)%vx2l,& - & dzero,p%wrk(level+1)%vx2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 + & vx2l,dzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') end if - - call inner_ml_aply(level+1,p,trans,work,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') + & a_err='Error during ADD smoother_apply') goto 9999 end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,p%wrk(level+1)%vy2l,& - & done,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + if (level < nlev) then + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + ! + ! Apply the prolongator + ! + 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,& + & a_err='Error during prolongation') + goto 9999 + end if - end if + end if + end associate + call psb_erractionrestore(err_act) return @@ -597,170 +603,174 @@ 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')) - if (level < nlev) then - ! - ! Apply the first smoother - ! The residual has been prepared before the recursive call. - ! + 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 + ! The residual has been prepared before the recursive call. + ! - if (pre) then - 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,& - & sweeps,work,info,init='Z') + if (pre) then + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & 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,& + & vx2l,dzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during PRE smoother_apply') + goto 9999 + end if + endif + ! + ! Compute the residual and call recursively + ! + if (pre) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + 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,& + & a_err='Error during restriction') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - endif - ! - ! 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) - - 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_) 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if - else - ! Shortcut: just transfer x2l. - call psb_map_X2Y(done,p%wrk(level)%vx2l,& - & dzero,p%wrk(level+1)%vx2l,& + ! Shortcut: just transfer x2l. + 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,& + & a_err='Error during restriction') + goto 9999 + end if + endif + + call inner_ml_aply(level+1,p,trans,work,info) + + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during restriction') + & a_err='Error during prolongation') goto 9999 end if - endif - call inner_ml_aply(level+1,p,trans,work,info) + if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then + + 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,& + & a_err='Error during W-cycle restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + + 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 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during W recusion/prolongation') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,p%wrk(level+1)%vy2l,& - & done,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + endif - 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W-cycle restriction') - goto 9999 - end if - - 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,& - & p%precv(level+1)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W recusion/prolongation') - goto 9999 - end if - - endif - - - 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,& - & 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 - ! - ! Apply the second smoother - ! - if (trans == 'N') then - sweeps = p%precv(level)%parms%sweeps_post - if (info == psb_success_) call p%precv(level)%sm2%apply(done,& - & p%wrk(level)%vty,done,p%wrk(level)%vy2l,& - & p%precv(level)%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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if - - endif - - else if (level == nlev) 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,& - & sweeps,work,info) + if (post) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + + ! + ! Apply the second smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & 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,& + & vty,done,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if - else + endif - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - end if + else if (level == nlev) then + + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & vx2l,dzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info) + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') + goto 9999 + end if + end associate + call psb_erractionrestore(err_act) return @@ -824,147 +834,150 @@ contains !K cycle - 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,& - & sweeps,work,info,init='Z') - - else if (level < nlev) then - - if (trans == 'N') then + 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,& - & 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 - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during 2-PRE smoother_apply') - goto 9999 - end if - + else if (level < nlev) then - ! - ! Compute the residual and call recursively - ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & 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,& + & vx2l,dzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if - 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_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during 2-PRE smoother_apply') + goto 9999 + end if - ! Apply the restriction - call psb_map_X2Y(done,p%wrk(level)%vty,& - & dzero,p%wrk(level + 1)%vx2l,& - & p%precv(level + 1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if + ! + ! Compute the residual and call recursively + ! - !Set the preconditioner + call psb_geaxpby(done,vx2l,& + & dzero,vty,& + & base_desc,info) - if (level <= nlev - 2 ) then - if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then - call mld_dinneritkcycle(p, level + 1, trans, work, 'FCG') - elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then - call mld_dinneritkcycle(p, level + 1, trans, work, 'GCR') - else + 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='Bad value for ml_cycle') + & a_err='Error during residue') goto 9999 + end if + + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + !Set the preconditioner + + if (level <= nlev - 2 ) then + if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then + call mld_dinneritkcycle(p, level + 1, trans, work, 'FCG') + elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then + call mld_dinneritkcycle(p, level + 1, trans, work, 'GCR') + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Bad value for ml_cycle') + goto 9999 + endif + else + call inner_ml_aply(level + 1 ,p,trans,work,info) endif - else - call inner_ml_aply(level + 1 ,p,trans,work,info) - endif - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,p%wrk(level+1)%vy2l,& - & done,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during prolongation') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - ! - ! 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,& - & 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 - ! - ! Apply the smoother - ! - 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,& - & sweeps,work,info,init='Z') + ! + ! Compute the residual + ! + 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,& + & a_err='Error during residue') + goto 9999 + end if + ! + ! Apply the smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & 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,& + & vty,done,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') goto 9999 - end if - else - - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - - endif + endif + end associate call psb_erractionrestore(err_act) return @@ -998,141 +1011,145 @@ contains real(psb_dpk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx character(len=20) :: name = 'innerit_k_cycle' - - !Assemble rhs, w, v, v1, x - - call psb_geasb(rhs,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(w,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v1,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(x,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - - if (psb_errstatus_fatal()) then - nc2l = p%precv(level)%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) - - !Apply the preconditioner - call p%wrk(level)%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_spmm(done,p%precv(level)%base_a,d(idx),dzero,v,p%precv(level)%base_desc,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if - - !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) - !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) - else - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid inner solver') - goto 9999 - endif - alpha = delta_old/tau - !Update residual w - call psb_geaxpby(-alpha, v, done, w, p%precv(level)%base_desc, info) + 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,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(w,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v1,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(x,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + call psb_geasb(d(1),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + + + call x%zero() + + ! rhs=vx2l and w=rhs + 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 = 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 - l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) - iter = 0 + delta0 = psb_genrm2(w, base_desc, info) - if (l2_norm <= rtol*delta0) then - !Update solution x - call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info) - else - iter = iter + 1 - idx=mod(iter,2) + !Apply the preconditioner + call vy2l%zero() - !Apply preconditioner - call psb_geaxpby(done,w,dzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info) + 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) - !Sparse matrix vector product + call psb_geaxpby(done,vy2l,dzero,d(idx),base_desc,info) - 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,v,base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during residue') goto 9999 end if - !tau1, tau2, tau3, tau4 + !FCG 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) - tau4= tau2 - (tau1*tau1)/tau + 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 - 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) - tau4= tau2 - (tau1*tau1)/tau + 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') goto 9999 endif - !Update solution - alpha=alpha-(tau1*tau3)/(tau*tau4) - call psb_geaxpby(alpha,d(idx - 1),done,x,p%precv(level)%base_desc,info) - alpha=tau3/tau4 - call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info) - endif + alpha = delta_old/tau + !Update residual w + call psb_geaxpby(-alpha, v, done, w, 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, base_desc, info) + else + iter = iter + 1 + idx=mod(iter,2) - call psb_geaxpby(done,x,dzero,p%wrk(level)%vy2l,p%precv(level)%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) + !Apply preconditioner + call psb_geaxpby(done,w,dzero,vx2l,base_desc,info) + call inner_ml_aply(level,p,trans,work,info) + call psb_geaxpby(done,vy2l,dzero,d(idx),base_desc,info) + + !Sparse matrix vector product + + 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') + goto 9999 + end if + + !tau1, tau2, tau3, tau4 + if (psb_toupper(trim(innersolv)) == 'FCG') then + 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, 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,& + & a_err='Invalid inner solver') + goto 9999 + endif + + !Update solution + alpha=alpha-(tau1*tau3)/(tau*tau4) + call psb_geaxpby(alpha,d(idx - 1),done,x,base_desc,info) + alpha=tau3/tau4 + call psb_geaxpby(alpha,d(idx),done,x,base_desc,info) + endif + call psb_geaxpby(done,x,dzero,vy2l,base_desc,info) + !Free vectors + 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 diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index 5b1ed296..7ab23307 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -251,20 +251,50 @@ 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 - ! - ! 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy') - goto 9999 - end if - - do isweep = 1, p%outer_sweeps - 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 + + + 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 + end if + + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call p%wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + 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,vx2l,base_desc,info) + call psb_spmm(-sone,base_a,y,sone,vx2l,base_desc,info) + end do + + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + ! ! With the current implementation, y2l is zeroed internally at first smoother. ! call p%wrk(level)%vy2l%zero() @@ -276,38 +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) - ! 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) - end do - - ! - ! If outer_sweeps == 1 we have just skipped the loop, and it's - ! equivalent to a single application. - ! - - ! - ! With the current implementation, y2l is zeroed internally at first smoother. - ! call p%wrk(level)%vy2l%zero() - ! - call inner_ml_aply(level,p,trans_,work,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,70 +480,75 @@ contains goto 9999 end if - 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) + 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,& + & 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,& + & vy2l,szero,vty,& + & base_desc, trans,& + & ione,work,info,init='Z') + + call p%precv(level)%sm2a%apply(sone,& + & vty,szero,vy2l,& + & base_desc, trans,& + & ione,work,info,init='Z') + end do - sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) - do k=1, sweeps + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(sone,& - & p%wrk(level)%vy2l,szero,p%wrk(level)%vtx,& - & p%precv(level)%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,& - & 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - ! Apply the restriction - call psb_map_X2Y(sone,p%wrk(level)%vx2l,& - & szero,p%wrk(level+1)%vx2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 + & vx2l,szero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') end if - - call inner_ml_aply(level+1,p,trans,work,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') + & a_err='Error during ADD smoother_apply') goto 9999 end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,& - & sone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + if (level < nlev) then + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + ! + ! Apply the prolongator + ! + 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,& + & a_err='Error during prolongation') + goto 9999 + end if - end if + end if + end associate + call psb_erractionrestore(err_act) return @@ -597,170 +603,174 @@ 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')) - if (level < nlev) then - ! - ! Apply the first smoother - ! The residual has been prepared before the recursive call. - ! + 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 + ! The residual has been prepared before the recursive call. + ! - if (pre) then - 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,& - & sweeps,work,info,init='Z') + if (pre) then + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(sone,& + & 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,& + & vx2l,szero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during PRE smoother_apply') + goto 9999 + end if + endif + ! + ! Compute the residual and call recursively + ! + if (pre) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + 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,& + & a_err='Error during restriction') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - endif - ! - ! 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) - - 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_) 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if - else - ! Shortcut: just transfer x2l. - call psb_map_X2Y(sone,p%wrk(level)%vx2l,& - & szero,p%wrk(level+1)%vx2l,& + ! Shortcut: just transfer x2l. + 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,& + & a_err='Error during restriction') + goto 9999 + end if + endif + + call inner_ml_aply(level+1,p,trans,work,info) + + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during restriction') + & a_err='Error during prolongation') goto 9999 end if - endif - call inner_ml_aply(level+1,p,trans,work,info) + if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then + + 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,& + & a_err='Error during W-cycle restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + + 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 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during W recusion/prolongation') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,& - & sone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + endif - 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W-cycle restriction') - goto 9999 - end if - - 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,& - & p%precv(level+1)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W recusion/prolongation') - goto 9999 - end if - - endif - - - 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,& - & 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 - ! - ! Apply the second smoother - ! - if (trans == 'N') then - sweeps = p%precv(level)%parms%sweeps_post - if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& - & p%wrk(level)%vty,sone,p%wrk(level)%vy2l,& - & p%precv(level)%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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if - - endif - - else if (level == nlev) 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,& - & sweeps,work,info) + if (post) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + + ! + ! Apply the second smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& + & 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,& + & vty,sone,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if - else + endif - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - end if + else if (level == nlev) then + + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(sone,& + & vx2l,szero,vy2l,& + & base_desc, trans,& + & sweeps,work,info) + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') + goto 9999 + end if + end associate + call psb_erractionrestore(err_act) return @@ -824,147 +834,150 @@ contains !K cycle - 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,& - & sweeps,work,info,init='Z') - - else if (level < nlev) then - - if (trans == 'N') then + 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,& - & 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 - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during 2-PRE smoother_apply') - goto 9999 - end if - + else if (level < nlev) then - ! - ! Compute the residual and call recursively - ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(sone,& + & 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,& + & vx2l,szero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if - 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_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during 2-PRE smoother_apply') + goto 9999 + end if - ! Apply the restriction - call psb_map_X2Y(sone,p%wrk(level)%vty,& - & szero,p%wrk(level + 1)%vx2l,& - & p%precv(level + 1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if + ! + ! Compute the residual and call recursively + ! - !Set the preconditioner + call psb_geaxpby(sone,vx2l,& + & szero,vty,& + & base_desc,info) - if (level <= nlev - 2 ) then - if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then - call mld_sinneritkcycle(p, level + 1, trans, work, 'FCG') - elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then - call mld_sinneritkcycle(p, level + 1, trans, work, 'GCR') - else + 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='Bad value for ml_cycle') + & a_err='Error during residue') goto 9999 + end if + + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + !Set the preconditioner + + if (level <= nlev - 2 ) then + if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then + call mld_sinneritkcycle(p, level + 1, trans, work, 'FCG') + elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then + call mld_sinneritkcycle(p, level + 1, trans, work, 'GCR') + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Bad value for ml_cycle') + goto 9999 + endif + else + call inner_ml_aply(level + 1 ,p,trans,work,info) endif - else - call inner_ml_aply(level + 1 ,p,trans,work,info) - endif - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,& - & sone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during prolongation') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - ! - ! 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,& - & 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 - ! - ! Apply the smoother - ! - 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,& - & sweeps,work,info,init='Z') + ! + ! Compute the residual + ! + 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,& + & a_err='Error during residue') + goto 9999 + end if + ! + ! Apply the smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& + & 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,& + & vty,sone,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') goto 9999 - end if - else - - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - - endif + endif + end associate call psb_erractionrestore(err_act) return @@ -998,141 +1011,145 @@ contains real(psb_spk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx character(len=20) :: name = 'innerit_k_cycle' - - !Assemble rhs, w, v, v1, x - - call psb_geasb(rhs,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(w,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v1,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(x,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - - if (psb_errstatus_fatal()) then - nc2l = p%precv(level)%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) - - !Apply the preconditioner - call p%wrk(level)%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_spmm(sone,p%precv(level)%base_a,d(idx),szero,v,p%precv(level)%base_desc,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if - - !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) - !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) - else - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid inner solver') - goto 9999 - endif - alpha = delta_old/tau - !Update residual w - call psb_geaxpby(-alpha, v, sone, w, p%precv(level)%base_desc, info) + 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,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(w,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v1,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(x,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + call psb_geasb(d(1),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + + + call x%zero() + + ! rhs=vx2l and w=rhs + 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 = 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 - l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) - iter = 0 + delta0 = psb_genrm2(w, base_desc, info) - if (l2_norm <= rtol*delta0) then - !Update solution x - call psb_geaxpby(alpha, d(idx), sone, x, p%precv(level)%base_desc, info) - else - iter = iter + 1 - idx=mod(iter,2) + !Apply the preconditioner + call vy2l%zero() - !Apply preconditioner - call psb_geaxpby(sone,w,szero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info) + 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) - !Sparse matrix vector product + call psb_geaxpby(sone,vy2l,szero,d(idx),base_desc,info) - 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,v,base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during residue') goto 9999 end if - !tau1, tau2, tau3, tau4 + !FCG 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) - tau4= tau2 - (tau1*tau1)/tau + 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 - 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) - tau4= tau2 - (tau1*tau1)/tau + 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') goto 9999 endif - !Update solution - alpha=alpha-(tau1*tau3)/(tau*tau4) - call psb_geaxpby(alpha,d(idx - 1),sone,x,p%precv(level)%base_desc,info) - alpha=tau3/tau4 - call psb_geaxpby(alpha,d(idx),sone,x,p%precv(level)%base_desc,info) - endif + alpha = delta_old/tau + !Update residual w + call psb_geaxpby(-alpha, v, sone, w, 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, base_desc, info) + else + iter = iter + 1 + idx=mod(iter,2) - call psb_geaxpby(sone,x,szero,p%wrk(level)%vy2l,p%precv(level)%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) + !Apply preconditioner + call psb_geaxpby(sone,w,szero,vx2l,base_desc,info) + call inner_ml_aply(level,p,trans,work,info) + call psb_geaxpby(sone,vy2l,szero,d(idx),base_desc,info) + + !Sparse matrix vector product + + 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') + goto 9999 + end if + + !tau1, tau2, tau3, tau4 + if (psb_toupper(trim(innersolv)) == 'FCG') then + 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, 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,& + & a_err='Invalid inner solver') + goto 9999 + endif + + !Update solution + alpha=alpha-(tau1*tau3)/(tau*tau4) + call psb_geaxpby(alpha,d(idx - 1),sone,x,base_desc,info) + alpha=tau3/tau4 + call psb_geaxpby(alpha,d(idx),sone,x,base_desc,info) + endif + call psb_geaxpby(sone,x,szero,vy2l,base_desc,info) + !Free vectors + 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 diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index bf7d7e48..372b6835 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -251,20 +251,50 @@ 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 - ! - ! 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy') - goto 9999 - end if - - do isweep = 1, p%outer_sweeps - 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 + + + 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 + end if + + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call p%wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + 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,vx2l,base_desc,info) + call psb_spmm(-zone,base_a,y,zone,vx2l,base_desc,info) + end do + + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + ! ! With the current implementation, y2l is zeroed internally at first smoother. ! call p%wrk(level)%vy2l%zero() @@ -276,38 +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) - ! 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) - end do - - ! - ! If outer_sweeps == 1 we have just skipped the loop, and it's - ! equivalent to a single application. - ! - - ! - ! With the current implementation, y2l is zeroed internally at first smoother. - ! call p%wrk(level)%vy2l%zero() - ! - call inner_ml_aply(level,p,trans_,work,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,70 +480,75 @@ contains goto 9999 end if - 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) + 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,& + & 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,& + & vy2l,zzero,vty,& + & base_desc, trans,& + & ione,work,info,init='Z') + + call p%precv(level)%sm2a%apply(zone,& + & vty,zzero,vy2l,& + & base_desc, trans,& + & ione,work,info,init='Z') + end do - sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) - do k=1, sweeps + else + sweeps = p%precv(level)%parms%sweeps_pre call p%precv(level)%sm%apply(zone,& - & p%wrk(level)%vy2l,zzero,p%wrk(level)%vtx,& - & p%precv(level)%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,& - & 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - ! Apply the restriction - call psb_map_X2Y(zone,p%wrk(level)%vx2l,& - & zzero,p%wrk(level+1)%vx2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 + & vx2l,zzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') end if - - call inner_ml_aply(level+1,p,trans,work,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') + & a_err='Error during ADD smoother_apply') goto 9999 end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,& - & zone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + if (level < nlev) then + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + ! + ! Apply the prolongator + ! + 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,& + & a_err='Error during prolongation') + goto 9999 + end if - end if + end if + end associate + call psb_erractionrestore(err_act) return @@ -597,170 +603,174 @@ 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')) - if (level < nlev) then - ! - ! Apply the first smoother - ! The residual has been prepared before the recursive call. - ! + 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 + ! The residual has been prepared before the recursive call. + ! - if (pre) then - 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,& - & sweeps,work,info,init='Z') + if (pre) then + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(zone,& + & 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,& + & vx2l,zzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during PRE smoother_apply') + goto 9999 + end if + endif + ! + ! Compute the residual and call recursively + ! + if (pre) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + 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,& + & a_err='Error during restriction') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - endif - ! - ! 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) - - 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_) 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if - else - ! Shortcut: just transfer x2l. - call psb_map_X2Y(zone,p%wrk(level)%vx2l,& - & zzero,p%wrk(level+1)%vx2l,& + ! Shortcut: just transfer x2l. + 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,& + & a_err='Error during restriction') + goto 9999 + end if + endif + + call inner_ml_aply(level+1,p,trans,work,info) + + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during restriction') + & a_err='Error during prolongation') goto 9999 end if - endif - call inner_ml_aply(level+1,p,trans,work,info) + if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then + + 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,& + & a_err='Error during W-cycle restriction') + goto 9999 + end if + + call inner_ml_aply(level+1,p,trans,work,info) + + 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 + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during W recusion/prolongation') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,& - & zone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + endif - 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,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W-cycle restriction') - goto 9999 - end if - - 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,& - & p%precv(level+1)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during W recusion/prolongation') - goto 9999 - end if - - endif - - - 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,& - & 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 - ! - ! Apply the second smoother - ! - if (trans == 'N') then - sweeps = p%precv(level)%parms%sweeps_post - if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& - & p%wrk(level)%vty,zone,p%wrk(level)%vy2l,& - & p%precv(level)%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,& - & sweeps,work,info,init='Z') - end if - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if - - endif - - else if (level == nlev) 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,& - & sweeps,work,info) + if (post) then + 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,& + & a_err='Error during residue') + goto 9999 + end if + + ! + ! Apply the second smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& + & 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,& + & vty,zone,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if - else + endif - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - end if + else if (level == nlev) then + + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(zone,& + & vx2l,zzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info) + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') + goto 9999 + end if + end associate + call psb_erractionrestore(err_act) return @@ -824,147 +834,150 @@ contains !K cycle - 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,& - & sweeps,work,info,init='Z') - - else if (level < nlev) then - - if (trans == 'N') then + 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,& - & 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 - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during 2-PRE smoother_apply') - goto 9999 - end if - + else if (level < nlev) then - ! - ! Compute the residual and call recursively - ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(zone,& + & 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,& + & vx2l,zzero,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if - 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_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during 2-PRE smoother_apply') + goto 9999 + end if - ! Apply the restriction - call psb_map_X2Y(zone,p%wrk(level)%vty,& - & zzero,p%wrk(level + 1)%vx2l,& - & p%precv(level + 1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if + ! + ! Compute the residual and call recursively + ! - !Set the preconditioner + call psb_geaxpby(zone,vx2l,& + & zzero,vty,& + & base_desc,info) - if (level <= nlev - 2 ) then - if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then - call mld_zinneritkcycle(p, level + 1, trans, work, 'FCG') - elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then - call mld_zinneritkcycle(p, level + 1, trans, work, 'GCR') - else + 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='Bad value for ml_cycle') + & a_err='Error during residue') goto 9999 + end if + + ! Apply the restriction + 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,& + & a_err='Error during restriction') + goto 9999 + end if + + !Set the preconditioner + + if (level <= nlev - 2 ) then + if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then + call mld_zinneritkcycle(p, level + 1, trans, work, 'FCG') + elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then + call mld_zinneritkcycle(p, level + 1, trans, work, 'GCR') + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Bad value for ml_cycle') + goto 9999 + endif + else + call inner_ml_aply(level + 1 ,p,trans,work,info) endif - else - call inner_ml_aply(level + 1 ,p,trans,work,info) - endif - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if - ! - ! Apply the prolongator - ! - call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,& - & zone,p%wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) + ! + ! Apply the prolongator + ! + 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,& - & a_err='Error during prolongation') - goto 9999 - end if + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - ! - ! 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,& - & 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 - ! - ! Apply the smoother - ! - 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,& - & sweeps,work,info,init='Z') + ! + ! Compute the residual + ! + 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,& + & a_err='Error during residue') + goto 9999 + end if + ! + ! Apply the smoother + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& + & 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,& + & vty,zone,vy2l,& + & base_desc, trans,& + & sweeps,work,info,init='Z') + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during POST smoother_apply') + goto 9999 + end if 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,& - & sweeps,work,info,init='Z') - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL vs NLEV') goto 9999 - end if - else - - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='Invalid LEVEL vs NLEV') - goto 9999 - - endif + endif + end associate call psb_erractionrestore(err_act) return @@ -998,141 +1011,145 @@ contains complex(psb_dpk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx character(len=20) :: name = 'innerit_k_cycle' - - !Assemble rhs, w, v, v1, x - - call psb_geasb(rhs,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(w,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(v1,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%vx2l%v) - call psb_geasb(x,& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=p%wrk(level)%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) - - if (psb_errstatus_fatal()) then - nc2l = p%precv(level)%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) - - !Apply the preconditioner - call p%wrk(level)%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_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v,p%precv(level)%base_desc,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if - - !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) - !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) - else - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid inner solver') - goto 9999 - endif - alpha = delta_old/tau - !Update residual w - call psb_geaxpby(-alpha, v, zone, w, p%precv(level)%base_desc, info) + 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,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(w,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(v1,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + call psb_geasb(x,& + & base_desc,info,& + & scratch=.true.,mold=vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + call psb_geasb(d(1),& + & base_desc,info,& + & scratch=.true.,mold=vy2l%v) + + + call x%zero() + + ! rhs=vx2l and w=rhs + 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 = 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 - l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) - iter = 0 + delta0 = psb_genrm2(w, base_desc, info) - if (l2_norm <= rtol*delta0) then - !Update solution x - call psb_geaxpby(alpha, d(idx), zone, x, p%precv(level)%base_desc, info) - else - iter = iter + 1 - idx=mod(iter,2) + !Apply the preconditioner + call vy2l%zero() - !Apply preconditioner - call psb_geaxpby(zone,w,zzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info) + 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) - !Sparse matrix vector product + call psb_geaxpby(zone,vy2l,zzero,d(idx),base_desc,info) - 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,v,base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during residue') goto 9999 end if - !tau1, tau2, tau3, tau4 + !FCG 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) - tau4= tau2 - (tau1*tau1)/tau + 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 - 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) - tau4= tau2 - (tau1*tau1)/tau + 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') goto 9999 endif - !Update solution - alpha=alpha-(tau1*tau3)/(tau*tau4) - call psb_geaxpby(alpha,d(idx - 1),zone,x,p%precv(level)%base_desc,info) - alpha=tau3/tau4 - call psb_geaxpby(alpha,d(idx),zone,x,p%precv(level)%base_desc,info) - endif + alpha = delta_old/tau + !Update residual w + call psb_geaxpby(-alpha, v, zone, w, 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, base_desc, info) + else + iter = iter + 1 + idx=mod(iter,2) - call psb_geaxpby(zone,x,zzero,p%wrk(level)%vy2l,p%precv(level)%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) + !Apply preconditioner + call psb_geaxpby(zone,w,zzero,vx2l,base_desc,info) + call inner_ml_aply(level,p,trans,work,info) + call psb_geaxpby(zone,vy2l,zzero,d(idx),base_desc,info) + + !Sparse matrix vector product + + 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') + goto 9999 + end if + + !tau1, tau2, tau3, tau4 + if (psb_toupper(trim(innersolv)) == 'FCG') then + 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, 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,& + & a_err='Invalid inner solver') + goto 9999 + endif + + !Update solution + alpha=alpha-(tau1*tau3)/(tau*tau4) + call psb_geaxpby(alpha,d(idx - 1),zone,x,base_desc,info) + alpha=tau3/tau4 + call psb_geaxpby(alpha,d(idx),zone,x,base_desc,info) + endif + call psb_geaxpby(zone,x,zzero,vy2l,base_desc,info) + !Free vectors + 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 diff --git a/mlprec/mld_c_onelev_mod.f90 b/mlprec/mld_c_onelev_mod.f90 index d42b2b0e..b76119a3 100644 --- a/mlprec/mld_c_onelev_mod.f90 +++ b/mlprec/mld_c_onelev_mod.f90 @@ -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 diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index bb6ae569..ce9d0157 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -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 diff --git a/mlprec/mld_s_onelev_mod.f90 b/mlprec/mld_s_onelev_mod.f90 index c81488ac..88f6bf14 100644 --- a/mlprec/mld_s_onelev_mod.f90 +++ b/mlprec/mld_s_onelev_mod.f90 @@ -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 diff --git a/mlprec/mld_z_onelev_mod.f90 b/mlprec/mld_z_onelev_mod.f90 index 619c4053..94a8fe43 100644 --- a/mlprec/mld_z_onelev_mod.f90 +++ b/mlprec/mld_z_onelev_mod.f90 @@ -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