diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index 618026e3..4825997f 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -703,6 +703,7 @@ contains if (level < nlev) then ! ! Apply the first smoother + ! The residual has been prepared before the recursive call. ! if (pre) then @@ -764,23 +765,8 @@ contains endif ! First guess is zero call mlprec_wrk(level+1)%vy2l%zero() - - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - - if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then - ! On second call will use output y2l as initial guess - if (info == psb_success_) & - & call inner_ml_aply(level+1,p,mlprec_wrk,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 - - + ! ! Apply the prolongator ! @@ -792,24 +778,54 @@ contains & a_err='Error during prolongation') goto 9999 end if + + if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then + ! On second call will use output y2l as initial guess + call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,& + & czero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + + if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info == psb_success_) call psb_map_X2Y(cone,mlprec_wrk(level)%vty,& + & czero,mlprec_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,mlprec_wrk,trans,work,info) + + if (info == psb_success_) call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& + & cone,mlprec_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 - if (.not.pre) then - ! - ! If we have only post, we need to compute the residual here. - ! - call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,& - & czero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) - call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& - & cone,mlprec_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_geaxpby(cone,mlprec_wrk(level)%vx2l,& + & czero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,& + & cone,mlprec_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 ! @@ -818,13 +834,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& & mlprec_wrk(level)%vty,cone,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(cone,& & mlprec_wrk(level)%vty,cone,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') end if if (info /= psb_success_) then diff --git a/mlprec/impl/mld_d_extprol_bld.f90 b/mlprec/impl/mld_d_extprol_bld.f90 index 27188e12..f36a6444 100644 --- a/mlprec/impl/mld_d_extprol_bld.f90 +++ b/mlprec/impl/mld_d_extprol_bld.f90 @@ -84,9 +84,9 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! Arguments type(psb_dspmat_type),intent(in), target :: a - type(psb_dspmat_type),intent(in), target :: prolv(:) - type(psb_dspmat_type),intent(in), target :: restrv(:) - type(psb_desc_type), intent(inout), target :: desc_a + type(psb_dspmat_type),intent(inout), target :: prolv(:) + type(psb_dspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a type(mld_dprec_type),intent(inout),target :: p integer(psb_ipk_), intent(out) :: info class(psb_d_base_sparse_mat), intent(in), optional :: amold @@ -303,7 +303,7 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) & mld_distr_mat_,is_distr_ml_coarse_mat) end if - call mld_d_bld_extaggr(p%precv(i-1)%base_a,& + call mld_d_extaggr_bld(p%precv(i-1)%base_a,& & p%precv(i-1)%base_desc,p%precv(i),restrv(i-1),prolv(i-1),info) diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index 1c48a8e3..a1a15480 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -528,7 +528,7 @@ contains call mld_d_inner_mult(p, mlprec_wrk, level, trans, work) case(mld_kcycle_ml_, mld_kcyclesym_ml_) - + call mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work) case default @@ -699,10 +699,11 @@ 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. ! if (pre) then @@ -719,7 +720,7 @@ contains & p%precv(level)%base_desc, trans,& & sweeps,work,info,init='Y') end if - + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during PRE smoother_apply') @@ -730,46 +731,42 @@ contains ! ! Compute the residual and call recursively ! - call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& - & dzero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) - - if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,done,mlprec_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,mlprec_wrk(level)%vty,& - & dzero,mlprec_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 - - + if (pre) then + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + + if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,done,mlprec_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,mlprec_wrk(level)%vty,& + & dzero,mlprec_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,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_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 + endif ! First guess is zero call mlprec_wrk(level+1)%vy2l%zero() - - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then - ! On second call will use output y2l as initial guess - if (info == psb_success_) & - & call inner_ml_aply(level+1,p,mlprec_wrk,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 - ! ! Apply the prolongator ! @@ -782,69 +779,80 @@ contains goto 9999 end if - + if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then + ! On second call will use output y2l as initial guess + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + + if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info == psb_success_) call psb_map_X2Y(done,mlprec_wrk(level)%vty,& + & dzero,mlprec_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,mlprec_wrk,trans,work,info) + + if (info == psb_success_) call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& + & done,mlprec_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 - if (.false.) then - call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& - & dzero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) - - if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,done,mlprec_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,& - & mlprec_wrk(level)%vty,done,mlprec_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,& - & mlprec_wrk(level)%vty,done,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Z') - end if - - else - - ! - ! Apply the second smoother - ! - if (trans == 'N') then - sweeps = p%precv(level)%parms%sweeps_post - if (info == psb_success_) call p%precv(level)%sm2%apply(done,& - & mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') - else - sweeps = p%precv(level)%parms%sweeps_pre - if (info == psb_success_) call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') - end if + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,& + & done,mlprec_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,& + & mlprec_wrk(level)%vty,done,mlprec_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,& + & mlprec_wrk(level)%vty,done,mlprec_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 - + + endif + else if (level == nlev) then - + sweeps = p%precv(level)%parms%sweeps if (info == psb_success_) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& @@ -892,8 +900,6 @@ contains integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: nlev, ilev, sweeps logical :: pre, post - real(psb_dpk_) :: tnrm - logical, parameter :: debug=.true. character(len=20) :: name @@ -956,10 +962,7 @@ contains & a_err='Error during 2-PRE smoother_apply') goto 9999 end if - if (debug) then - tnrm = psb_genrm2(mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) - write(0,*)' After first smoother y2l ',tnrm - end if + ! ! Compute the residual and call recursively @@ -977,20 +980,11 @@ contains & a_err='Error during residue') goto 9999 end if - if (debug) then - tnrm = psb_genrm2(mlprec_wrk(level)%vty,p%precv(level)%base_desc,info) - write(0,*)' Residual before restriction ',tnrm - end if ! Apply the restriction call psb_map_X2Y(done,mlprec_wrk(level)%vty,& & dzero,mlprec_wrk(level + 1)%vx2l,& & p%precv(level + 1)%map,info,work=work) - if (debug) then - tnrm = psb_genrm2(mlprec_wrk(level+1)%vx2l,p%precv(level+1)%base_desc,info) - write(0,*)' Output of restriction ',tnrm - end if - call mlprec_wrk(level + 1)%vy2l%zero() if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -1000,7 +994,7 @@ contains !Set the preconditioner - if ((level < nlev - 1)) then + if ((level < nlev - 2)) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then @@ -1030,21 +1024,22 @@ contains end if ! - ! Apply the smoother + ! Compute the residual ! call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& & dzero,mlprec_wrk(level)%vty,& & p%precv(level)%base_desc,info) - - if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info,work=work,trans=trans) + call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& + & done,mlprec_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 - + 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,& @@ -1105,7 +1100,6 @@ contains real(psb_dpk_) :: l2_norm, delta, rtol=0.25 real(psb_dpk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx - logical, parameter :: debug=.true. !Assemble rhs, w, v, v1, x @@ -1125,7 +1119,7 @@ contains & p%precv(level)%base_desc,info,& & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) -!!$ call x%zero() + call x%zero() ! rhs=vx2l and w=rhs call psb_geaxpby(done,mlprec_wrk(level)%vx2l,dzero,rhs,& @@ -1141,8 +1135,8 @@ contains goto 9999 end if - delta = psb_genrm2(w, p%precv(level)%base_desc, info) - if (debug) write(0,*)'On entry delta ',delta + delta = psb_gedot(w, w, p%precv(level)%base_desc, info) + !Apply the preconditioner call mlprec_wrk(level)%vy2l%set(dzero) @@ -1185,11 +1179,10 @@ contains l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info) iter = 0 - if (debug) write(0,*) 'Alpha ',alpha,' l2_norm',l2_norm,& - & ' delta_old',delta_old,' tau',tau + if (l2_norm <= rtol*delta) then !Update solution x - call psb_geaxpby(alpha, d(idx), dzero, x, p%precv(level)%base_desc, info) + call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info) else iter = iter + 1 idx=mod(iter,2) @@ -1225,7 +1218,7 @@ contains !Update solution alpha=alpha-(tau1*tau3)/(tau*tau4) - call psb_geaxpby(alpha,d(idx - 1),dzero,x,p%precv(level)%base_desc,info) + 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 diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index 43ebaf34..98afa860 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -703,6 +703,7 @@ contains if (level < nlev) then ! ! Apply the first smoother + ! The residual has been prepared before the recursive call. ! if (pre) then @@ -764,23 +765,8 @@ contains endif ! First guess is zero call mlprec_wrk(level+1)%vy2l%zero() - - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - - if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then - ! On second call will use output y2l as initial guess - if (info == psb_success_) & - & call inner_ml_aply(level+1,p,mlprec_wrk,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 - - + ! ! Apply the prolongator ! @@ -792,24 +778,54 @@ contains & a_err='Error during prolongation') goto 9999 end if + + if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then + ! On second call will use output y2l as initial guess + call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,& + & szero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + + if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info == psb_success_) call psb_map_X2Y(sone,mlprec_wrk(level)%vty,& + & szero,mlprec_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,mlprec_wrk,trans,work,info) + + if (info == psb_success_) call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& + & sone,mlprec_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 - if (.not.pre) then - ! - ! If we have only post, we need to compute the residual here. - ! - call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,& - & szero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) - call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& - & sone,mlprec_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_geaxpby(sone,mlprec_wrk(level)%vx2l,& + & szero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,& + & sone,mlprec_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 ! @@ -818,13 +834,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& & mlprec_wrk(level)%vty,sone,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(sone,& & mlprec_wrk(level)%vty,sone,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') end if if (info /= psb_success_) then diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index e5832b8e..50c50d68 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -703,6 +703,7 @@ contains if (level < nlev) then ! ! Apply the first smoother + ! The residual has been prepared before the recursive call. ! if (pre) then @@ -764,23 +765,8 @@ contains endif ! First guess is zero call mlprec_wrk(level+1)%vy2l%zero() - - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - - if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then - ! On second call will use output y2l as initial guess - if (info == psb_success_) & - & call inner_ml_aply(level+1,p,mlprec_wrk,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 - - + ! ! Apply the prolongator ! @@ -792,24 +778,54 @@ contains & a_err='Error during prolongation') goto 9999 end if + + if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then + ! On second call will use output y2l as initial guess + call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,& + & zzero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + + if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info == psb_success_) call psb_map_X2Y(zone,mlprec_wrk(level)%vty,& + & zzero,mlprec_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,mlprec_wrk,trans,work,info) + + if (info == psb_success_) call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& + & zone,mlprec_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 - if (.not.pre) then - ! - ! If we have only post, we need to compute the residual here. - ! - call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,& - & zzero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) - call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& - & zone,mlprec_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_geaxpby(zone,mlprec_wrk(level)%vx2l,& + & zzero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,& + & zone,mlprec_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 ! @@ -818,13 +834,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& & mlprec_wrk(level)%vty,zone,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(zone,& & mlprec_wrk(level)%vty,zone,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') end if if (info /= psb_success_) then diff --git a/mlprec/mld_d_prec_mod.f90 b/mlprec/mld_d_prec_mod.f90 index f90e974b..48919fc3 100644 --- a/mlprec/mld_d_prec_mod.f90 +++ b/mlprec/mld_d_prec_mod.f90 @@ -112,9 +112,9 @@ module mld_d_prec_mod implicit none ! Arguments type(psb_dspmat_type),intent(in), target :: a - type(psb_dspmat_type),intent(in), target :: prolv(:) - type(psb_dspmat_type),intent(in), target :: restrv(:) - type(psb_desc_type), intent(inout), target :: desc_a + type(psb_dspmat_type),intent(inout), target :: prolv(:) + type(psb_dspmat_type),intent(inout), target :: restrv(:) + type(psb_desc_type), intent(inout), target :: desc_a type(mld_dprec_type),intent(inout),target :: p integer(psb_ipk_), intent(out) :: info class(psb_d_base_sparse_mat), intent(in), optional :: amold