From 823db4f943aa7be9ef193195b76bfe070f6bb551 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 11 Dec 2017 10:21:00 +0000 Subject: [PATCH] Fixed allocte_wrk & free_wrk for WV allocation. Modified interface of smoothers to use WV. Initial tests. Added WV to calls to MAP_X2Y & MAP_Y2X. --- mlprec/impl/mld_cmlprec_aply.f90 | 68 ++++--- mlprec/impl/mld_cprecaply.f90 | 155 ++++++++-------- mlprec/impl/mld_cprecset.F90 | 4 +- mlprec/impl/mld_dmlprec_aply.f90 | 68 ++++--- mlprec/impl/mld_dprecaply.f90 | 155 ++++++++-------- mlprec/impl/mld_dprecset.F90 | 4 +- mlprec/impl/mld_smlprec_aply.f90 | 68 ++++--- mlprec/impl/mld_sprecaply.f90 | 155 ++++++++-------- mlprec/impl/mld_sprecset.F90 | 4 +- mlprec/impl/mld_zmlprec_aply.f90 | 68 ++++--- mlprec/impl/mld_zprecaply.f90 | 155 ++++++++-------- mlprec/impl/mld_zprecset.F90 | 4 +- .../smoother/mld_c_as_smoother_apply_vect.f90 | 167 +++++++++--------- .../mld_c_base_smoother_apply_vect.f90 | 4 +- .../mld_c_jac_smoother_apply_vect.f90 | 133 +++++++------- .../smoother/mld_d_as_smoother_apply_vect.f90 | 167 +++++++++--------- .../mld_d_base_smoother_apply_vect.f90 | 4 +- .../mld_d_jac_smoother_apply_vect.f90 | 133 +++++++------- .../smoother/mld_s_as_smoother_apply_vect.f90 | 167 +++++++++--------- .../mld_s_base_smoother_apply_vect.f90 | 4 +- .../mld_s_jac_smoother_apply_vect.f90 | 133 +++++++------- .../smoother/mld_z_as_smoother_apply_vect.f90 | 167 +++++++++--------- .../mld_z_base_smoother_apply_vect.f90 | 4 +- .../mld_z_jac_smoother_apply_vect.f90 | 133 +++++++------- mlprec/mld_c_as_smoother.f90 | 4 +- mlprec/mld_c_base_smoother_mod.f90 | 4 +- mlprec/mld_c_jac_smoother.f90 | 4 +- mlprec/mld_c_onelev_mod.f90 | 17 +- mlprec/mld_d_as_smoother.f90 | 4 +- mlprec/mld_d_base_smoother_mod.f90 | 4 +- mlprec/mld_d_jac_smoother.f90 | 4 +- mlprec/mld_d_onelev_mod.f90 | 17 +- mlprec/mld_s_as_smoother.f90 | 4 +- mlprec/mld_s_base_smoother_mod.f90 | 4 +- mlprec/mld_s_jac_smoother.f90 | 4 +- mlprec/mld_s_onelev_mod.f90 | 17 +- mlprec/mld_z_as_smoother.f90 | 4 +- mlprec/mld_z_base_smoother_mod.f90 | 4 +- mlprec/mld_z_jac_smoother.f90 | 4 +- mlprec/mld_z_onelev_mod.f90 | 17 +- 40 files changed, 1212 insertions(+), 1028 deletions(-) diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index 11629f69..3ba3d8c4 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -255,7 +255,8 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv(:)) ! ! At first iteration we must use the input BETA ! @@ -482,7 +483,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (allocated(p%precv(level)%sm2a)) then call psb_geaxpby(cone,& @@ -494,12 +496,12 @@ contains call p%precv(level)%sm%apply(cone,& & vy2l,czero,vty,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') call p%precv(level)%sm2a%apply(cone,& & vty,czero,vy2l,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') end do else @@ -507,7 +509,7 @@ contains call p%precv(level)%sm%apply(cone,& & vx2l,czero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -519,7 +521,8 @@ contains ! Apply the restriction call psb_map_X2Y(cone,vx2l,& & czero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -538,7 +541,8 @@ contains ! call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,& & cone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -605,7 +609,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level < nlev) then ! ! Apply the first smoother @@ -618,13 +623,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(cone,& & vx2l,czero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -651,7 +656,8 @@ contains end if call psb_map_X2Y(cone,vty,& & czero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -661,7 +667,8 @@ contains ! 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -676,7 +683,8 @@ contains ! call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,& & cone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -693,7 +701,8 @@ contains & 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -704,7 +713,8 @@ contains 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -737,13 +747,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& & vty,cone,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -760,7 +770,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(cone,& & vx2l,czero,vy2l,& & base_desc, trans,& - & sweeps,work,info) + & sweeps,work,wv,info) else @@ -836,7 +846,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level == nlev) then ! ! Apply smoother @@ -845,7 +856,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(cone,& & vx2l,czero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') else if (level < nlev) then @@ -854,13 +865,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(cone,& & vx2l,czero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -890,7 +901,8 @@ contains ! Apply the restriction call psb_map_X2Y(cone,vty,& & czero,p%precv(level + 1)%wrk%vx2l,& - & p%precv(level + 1)%map,info,work=work) + & p%precv(level + 1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -925,7 +937,8 @@ contains ! call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,& & cone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -955,13 +968,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& & vty,cone,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -1014,7 +1027,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) !Assemble rhs, w, v, v1, x diff --git a/mlprec/impl/mld_cprecaply.f90 b/mlprec/impl/mld_cprecaply.f90 index 4e8a99c4..3641d26e 100644 --- a/mlprec/impl/mld_cprecaply.f90 +++ b/mlprec/impl/mld_cprecaply.f90 @@ -323,6 +323,7 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work) complex(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_cprecaply' @@ -358,6 +359,10 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + if (size(prec%precv) >1) then ! ! Number of levels > 1: apply the multilevel preconditioner @@ -375,31 +380,29 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then - ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. - ! - twoside: block - type(psb_c_vect_type) :: w1,w2 - call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.) + associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,& + & wv => prec%precv(1)%wrk%wv) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! call psb_geaxpby(cone,x,czero,w1,desc_data,info) select case(trans_) case ('N') do k=1, nswps call prec%precv(1)%sm%apply(cone,w1,czero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm2a%apply(cone,w2,czero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do - + case('T','C') do k=1, nswps call prec%precv(1)%sm2a%apply(cone,w1,czero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm%apply(cone,w2,czero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do case default info = psb_err_from_subroutine_ @@ -407,14 +410,12 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work) goto 9999 end select call psb_geaxpby(cone,w1,czero,y,desc_data,info) - call psb_gefree(w1,desc_data,info) - call psb_gefree(w2,desc_data,info) - end block twoside - - else - call prec%precv(1)%sm%apply(cone,x,czero,y,desc_data,trans_,& - & nswps,work_,info) - end if + else + call prec%precv(1)%sm%apply(cone,x,czero,y,desc_data,trans_,& + & nswps,work_,wv,info) + end if + end associate + else info = psb_err_from_subroutine_ai_ @@ -426,6 +427,7 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! If the original distribution has an overlap we should fix that. call psb_halo(y,desc_data,info,data=psb_comm_mov_) + if (do_alloc_wrk) call prec%free_wrk(info) if (present(work)) then else @@ -459,10 +461,10 @@ subroutine mld_cprecaply1_vect(prec,x,desc_data,info,trans,work) ! Local variables character :: trans_ - type(psb_c_vect_type) :: ww complex(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_cprecaply' @@ -498,63 +500,68 @@ subroutine mld_cprecaply1_vect(prec,x,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if - call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) - if (size(prec%precv) >1) then - ! - ! Number of levels > 1: apply the multilevel preconditioner - ! - call mld_mlprec_aply(cone,prec,x,czero,ww,desc_data,trans_,work_,info) - if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_cmlprec_aply') - goto 9999 - end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) - else if (size(prec%precv) == 1) then - ! - ! Number of levels = 1: apply the base preconditioner - ! - nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then + associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv) + call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) + if (size(prec%precv) >1) then ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. + ! Number of levels > 1: apply the multilevel preconditioner ! - select case(trans_) - case ('N') - do k=1, nswps - call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm2a%apply(cone,ww,czero,x,desc_data,trans_,& - & ione, work_,info) - end do - case('T','C') - do k=1, nswps - call prec%precv(1)%sm2a%apply(cone,x,czero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm%apply(cone,ww,czero,x,desc_data,trans_,& - & ione, work_,info) - end do - case default - info = psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='Invalid trans') - goto 9999 - end select - - else - call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,& - & nswps, work_,info) + call mld_mlprec_aply(cone,prec,x,czero,ww,desc_data,trans_,work_,info) if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info) - end if - else + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_cmlprec_aply') + goto 9999 + end if - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='Invalid size of precv',& - & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) - goto 9999 - endif + else if (size(prec%precv) == 1) then + ! + ! Number of levels = 1: apply the base preconditioner + ! + nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! + select case(trans_) + case ('N') + do k=1, nswps + call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm2a%apply(cone,ww,czero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case('T','C') + do k=1, nswps + call prec%precv(1)%sm2a%apply(cone,x,czero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm%apply(cone,ww,czero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case default + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Invalid trans') + goto 9999 + end select - if (info == 0) call psb_gefree(ww,desc_data,info) + else + call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,& + & nswps, work_,wv,info) + if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info) + end if + else + + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='Invalid size of precv',& + & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) + goto 9999 + endif + end associate +!!$ if (info == 0) call psb_gefree(ww,desc_data,info) ! If the original distribution has an overlap we should fix that. call psb_halo(x,desc_data,info,data=psb_comm_mov_) diff --git a/mlprec/impl/mld_cprecset.F90 b/mlprec/impl/mld_cprecset.F90 index 0ef1fb91..307c4e81 100644 --- a/mlprec/impl/mld_cprecset.F90 +++ b/mlprec/impl/mld_cprecset.F90 @@ -432,7 +432,7 @@ subroutine mld_cprecsetsm(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsm' info = psb_success_ @@ -495,7 +495,7 @@ subroutine mld_cprecsetsv(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsv' info = psb_success_ diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index 4a046d80..e7d74c11 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -255,7 +255,8 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv(:)) ! ! At first iteration we must use the input BETA ! @@ -482,7 +483,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (allocated(p%precv(level)%sm2a)) then call psb_geaxpby(done,& @@ -494,12 +496,12 @@ contains call p%precv(level)%sm%apply(done,& & vy2l,dzero,vty,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') call p%precv(level)%sm2a%apply(done,& & vty,dzero,vy2l,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') end do else @@ -507,7 +509,7 @@ contains call p%precv(level)%sm%apply(done,& & vx2l,dzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -519,7 +521,8 @@ contains ! Apply the restriction call psb_map_X2Y(done,vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -538,7 +541,8 @@ contains ! call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,& & done,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -605,7 +609,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level < nlev) then ! ! Apply the first smoother @@ -618,13 +623,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(done,& & vx2l,dzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -651,7 +656,8 @@ contains end if call psb_map_X2Y(done,vty,& & dzero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -661,7 +667,8 @@ contains ! 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -676,7 +683,8 @@ contains ! call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,& & done,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -693,7 +701,8 @@ contains & 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -704,7 +713,8 @@ contains 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -737,13 +747,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(done,& & vty,done,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -760,7 +770,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(done,& & vx2l,dzero,vy2l,& & base_desc, trans,& - & sweeps,work,info) + & sweeps,work,wv,info) else @@ -836,7 +846,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level == nlev) then ! ! Apply smoother @@ -845,7 +856,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(done,& & vx2l,dzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') else if (level < nlev) then @@ -854,13 +865,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(done,& & vx2l,dzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -890,7 +901,8 @@ contains ! Apply the restriction call psb_map_X2Y(done,vty,& & dzero,p%precv(level + 1)%wrk%vx2l,& - & p%precv(level + 1)%map,info,work=work) + & p%precv(level + 1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -925,7 +937,8 @@ contains ! call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,& & done,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -955,13 +968,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(done,& & vty,done,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -1014,7 +1027,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) !Assemble rhs, w, v, v1, x diff --git a/mlprec/impl/mld_dprecaply.f90 b/mlprec/impl/mld_dprecaply.f90 index bc91bb66..d9277daf 100644 --- a/mlprec/impl/mld_dprecaply.f90 +++ b/mlprec/impl/mld_dprecaply.f90 @@ -323,6 +323,7 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work) real(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_dprecaply' @@ -358,6 +359,10 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + if (size(prec%precv) >1) then ! ! Number of levels > 1: apply the multilevel preconditioner @@ -375,31 +380,29 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then - ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. - ! - twoside: block - type(psb_d_vect_type) :: w1,w2 - call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.) + associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,& + & wv => prec%precv(1)%wrk%wv) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! call psb_geaxpby(done,x,dzero,w1,desc_data,info) select case(trans_) case ('N') do k=1, nswps call prec%precv(1)%sm%apply(done,w1,dzero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm2a%apply(done,w2,dzero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do - + case('T','C') do k=1, nswps call prec%precv(1)%sm2a%apply(done,w1,dzero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm%apply(done,w2,dzero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do case default info = psb_err_from_subroutine_ @@ -407,14 +410,12 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work) goto 9999 end select call psb_geaxpby(done,w1,dzero,y,desc_data,info) - call psb_gefree(w1,desc_data,info) - call psb_gefree(w2,desc_data,info) - end block twoside - - else - call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_,& - & nswps,work_,info) - end if + else + call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_,& + & nswps,work_,wv,info) + end if + end associate + else info = psb_err_from_subroutine_ai_ @@ -426,6 +427,7 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! If the original distribution has an overlap we should fix that. call psb_halo(y,desc_data,info,data=psb_comm_mov_) + if (do_alloc_wrk) call prec%free_wrk(info) if (present(work)) then else @@ -459,10 +461,10 @@ subroutine mld_dprecaply1_vect(prec,x,desc_data,info,trans,work) ! Local variables character :: trans_ - type(psb_d_vect_type) :: ww real(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_dprecaply' @@ -498,63 +500,68 @@ subroutine mld_dprecaply1_vect(prec,x,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if - call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) - if (size(prec%precv) >1) then - ! - ! Number of levels > 1: apply the multilevel preconditioner - ! - call mld_mlprec_aply(done,prec,x,dzero,ww,desc_data,trans_,work_,info) - if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_dmlprec_aply') - goto 9999 - end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) - else if (size(prec%precv) == 1) then - ! - ! Number of levels = 1: apply the base preconditioner - ! - nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then + associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv) + call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) + if (size(prec%precv) >1) then ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. + ! Number of levels > 1: apply the multilevel preconditioner ! - select case(trans_) - case ('N') - do k=1, nswps - call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm2a%apply(done,ww,dzero,x,desc_data,trans_,& - & ione, work_,info) - end do - case('T','C') - do k=1, nswps - call prec%precv(1)%sm2a%apply(done,x,dzero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm%apply(done,ww,dzero,x,desc_data,trans_,& - & ione, work_,info) - end do - case default - info = psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='Invalid trans') - goto 9999 - end select - - else - call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,& - & nswps, work_,info) + call mld_mlprec_aply(done,prec,x,dzero,ww,desc_data,trans_,work_,info) if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info) - end if - else + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_dmlprec_aply') + goto 9999 + end if - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='Invalid size of precv',& - & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) - goto 9999 - endif + else if (size(prec%precv) == 1) then + ! + ! Number of levels = 1: apply the base preconditioner + ! + nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! + select case(trans_) + case ('N') + do k=1, nswps + call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm2a%apply(done,ww,dzero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case('T','C') + do k=1, nswps + call prec%precv(1)%sm2a%apply(done,x,dzero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm%apply(done,ww,dzero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case default + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Invalid trans') + goto 9999 + end select - if (info == 0) call psb_gefree(ww,desc_data,info) + else + call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,& + & nswps, work_,wv,info) + if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info) + end if + else + + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='Invalid size of precv',& + & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) + goto 9999 + endif + end associate +!!$ if (info == 0) call psb_gefree(ww,desc_data,info) ! If the original distribution has an overlap we should fix that. call psb_halo(x,desc_data,info,data=psb_comm_mov_) diff --git a/mlprec/impl/mld_dprecset.F90 b/mlprec/impl/mld_dprecset.F90 index 3c74b628..9bf541e0 100644 --- a/mlprec/impl/mld_dprecset.F90 +++ b/mlprec/impl/mld_dprecset.F90 @@ -465,7 +465,7 @@ subroutine mld_dprecsetsm(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsm' info = psb_success_ @@ -528,7 +528,7 @@ subroutine mld_dprecsetsv(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsv' info = psb_success_ diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index 7ab23307..5a62c214 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -255,7 +255,8 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv(:)) ! ! At first iteration we must use the input BETA ! @@ -482,7 +483,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (allocated(p%precv(level)%sm2a)) then call psb_geaxpby(sone,& @@ -494,12 +496,12 @@ contains call p%precv(level)%sm%apply(sone,& & vy2l,szero,vty,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') call p%precv(level)%sm2a%apply(sone,& & vty,szero,vy2l,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') end do else @@ -507,7 +509,7 @@ contains call p%precv(level)%sm%apply(sone,& & vx2l,szero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -519,7 +521,8 @@ contains ! Apply the restriction call psb_map_X2Y(sone,vx2l,& & szero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -538,7 +541,8 @@ contains ! call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,& & sone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -605,7 +609,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level < nlev) then ! ! Apply the first smoother @@ -618,13 +623,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(sone,& & vx2l,szero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -651,7 +656,8 @@ contains end if call psb_map_X2Y(sone,vty,& & szero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -661,7 +667,8 @@ contains ! 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -676,7 +683,8 @@ contains ! call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,& & sone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -693,7 +701,8 @@ contains & 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -704,7 +713,8 @@ contains 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -737,13 +747,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& & vty,sone,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -760,7 +770,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(sone,& & vx2l,szero,vy2l,& & base_desc, trans,& - & sweeps,work,info) + & sweeps,work,wv,info) else @@ -836,7 +846,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level == nlev) then ! ! Apply smoother @@ -845,7 +856,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(sone,& & vx2l,szero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') else if (level < nlev) then @@ -854,13 +865,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(sone,& & vx2l,szero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -890,7 +901,8 @@ contains ! Apply the restriction call psb_map_X2Y(sone,vty,& & szero,p%precv(level + 1)%wrk%vx2l,& - & p%precv(level + 1)%map,info,work=work) + & p%precv(level + 1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -925,7 +937,8 @@ contains ! call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,& & sone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -955,13 +968,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& & vty,sone,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -1014,7 +1027,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) !Assemble rhs, w, v, v1, x diff --git a/mlprec/impl/mld_sprecaply.f90 b/mlprec/impl/mld_sprecaply.f90 index b68d3421..e4065b7b 100644 --- a/mlprec/impl/mld_sprecaply.f90 +++ b/mlprec/impl/mld_sprecaply.f90 @@ -323,6 +323,7 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work) real(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_sprecaply' @@ -358,6 +359,10 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + if (size(prec%precv) >1) then ! ! Number of levels > 1: apply the multilevel preconditioner @@ -375,31 +380,29 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then - ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. - ! - twoside: block - type(psb_s_vect_type) :: w1,w2 - call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.) + associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,& + & wv => prec%precv(1)%wrk%wv) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! call psb_geaxpby(sone,x,szero,w1,desc_data,info) select case(trans_) case ('N') do k=1, nswps call prec%precv(1)%sm%apply(sone,w1,szero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm2a%apply(sone,w2,szero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do - + case('T','C') do k=1, nswps call prec%precv(1)%sm2a%apply(sone,w1,szero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm%apply(sone,w2,szero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do case default info = psb_err_from_subroutine_ @@ -407,14 +410,12 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work) goto 9999 end select call psb_geaxpby(sone,w1,szero,y,desc_data,info) - call psb_gefree(w1,desc_data,info) - call psb_gefree(w2,desc_data,info) - end block twoside - - else - call prec%precv(1)%sm%apply(sone,x,szero,y,desc_data,trans_,& - & nswps,work_,info) - end if + else + call prec%precv(1)%sm%apply(sone,x,szero,y,desc_data,trans_,& + & nswps,work_,wv,info) + end if + end associate + else info = psb_err_from_subroutine_ai_ @@ -426,6 +427,7 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! If the original distribution has an overlap we should fix that. call psb_halo(y,desc_data,info,data=psb_comm_mov_) + if (do_alloc_wrk) call prec%free_wrk(info) if (present(work)) then else @@ -459,10 +461,10 @@ subroutine mld_sprecaply1_vect(prec,x,desc_data,info,trans,work) ! Local variables character :: trans_ - type(psb_s_vect_type) :: ww real(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_sprecaply' @@ -498,63 +500,68 @@ subroutine mld_sprecaply1_vect(prec,x,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if - call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) - if (size(prec%precv) >1) then - ! - ! Number of levels > 1: apply the multilevel preconditioner - ! - call mld_mlprec_aply(sone,prec,x,szero,ww,desc_data,trans_,work_,info) - if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_smlprec_aply') - goto 9999 - end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) - else if (size(prec%precv) == 1) then - ! - ! Number of levels = 1: apply the base preconditioner - ! - nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then + associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv) + call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) + if (size(prec%precv) >1) then ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. + ! Number of levels > 1: apply the multilevel preconditioner ! - select case(trans_) - case ('N') - do k=1, nswps - call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm2a%apply(sone,ww,szero,x,desc_data,trans_,& - & ione, work_,info) - end do - case('T','C') - do k=1, nswps - call prec%precv(1)%sm2a%apply(sone,x,szero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm%apply(sone,ww,szero,x,desc_data,trans_,& - & ione, work_,info) - end do - case default - info = psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='Invalid trans') - goto 9999 - end select - - else - call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,& - & nswps, work_,info) + call mld_mlprec_aply(sone,prec,x,szero,ww,desc_data,trans_,work_,info) if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info) - end if - else + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_smlprec_aply') + goto 9999 + end if - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='Invalid size of precv',& - & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) - goto 9999 - endif + else if (size(prec%precv) == 1) then + ! + ! Number of levels = 1: apply the base preconditioner + ! + nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! + select case(trans_) + case ('N') + do k=1, nswps + call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm2a%apply(sone,ww,szero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case('T','C') + do k=1, nswps + call prec%precv(1)%sm2a%apply(sone,x,szero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm%apply(sone,ww,szero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case default + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Invalid trans') + goto 9999 + end select - if (info == 0) call psb_gefree(ww,desc_data,info) + else + call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,& + & nswps, work_,wv,info) + if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info) + end if + else + + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='Invalid size of precv',& + & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) + goto 9999 + endif + end associate +!!$ if (info == 0) call psb_gefree(ww,desc_data,info) ! If the original distribution has an overlap we should fix that. call psb_halo(x,desc_data,info,data=psb_comm_mov_) diff --git a/mlprec/impl/mld_sprecset.F90 b/mlprec/impl/mld_sprecset.F90 index 688757d4..f8507f94 100644 --- a/mlprec/impl/mld_sprecset.F90 +++ b/mlprec/impl/mld_sprecset.F90 @@ -432,7 +432,7 @@ subroutine mld_sprecsetsm(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsm' info = psb_success_ @@ -495,7 +495,7 @@ subroutine mld_sprecsetsv(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsv' info = psb_success_ diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index 372b6835..920ba17f 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -255,7 +255,8 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv(:)) ! ! At first iteration we must use the input BETA ! @@ -482,7 +483,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (allocated(p%precv(level)%sm2a)) then call psb_geaxpby(zone,& @@ -494,12 +496,12 @@ contains call p%precv(level)%sm%apply(zone,& & vy2l,zzero,vty,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') call p%precv(level)%sm2a%apply(zone,& & vty,zzero,vy2l,& & base_desc, trans,& - & ione,work,info,init='Z') + & ione,work,wv,info,init='Z') end do else @@ -507,7 +509,7 @@ contains call p%precv(level)%sm%apply(zone,& & vx2l,zzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -519,7 +521,8 @@ contains ! Apply the restriction call psb_map_X2Y(zone,vx2l,& & zzero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -538,7 +541,8 @@ contains ! call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,& & zone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -605,7 +609,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level < nlev) then ! ! Apply the first smoother @@ -618,13 +623,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(zone,& & vx2l,zzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -651,7 +656,8 @@ contains end if call psb_map_X2Y(zone,vty,& & zzero,p%precv(level+1)%wrk%vx2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -661,7 +667,8 @@ contains ! 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -676,7 +683,8 @@ contains ! call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,& & zone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -693,7 +701,8 @@ contains & 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -704,7 +713,8 @@ contains 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) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -737,13 +747,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& & vty,zone,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -760,7 +770,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(zone,& & vx2l,zzero,vy2l,& & base_desc, trans,& - & sweeps,work,info) + & sweeps,work,wv,info) else @@ -836,7 +846,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) if (level == nlev) then ! ! Apply smoother @@ -845,7 +856,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(zone,& & vx2l,zzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,info,init='Z') else if (level < nlev) then @@ -854,13 +865,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(zone,& & vx2l,zzero,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -890,7 +901,8 @@ contains ! Apply the restriction call psb_map_X2Y(zone,vty,& & zzero,p%precv(level + 1)%wrk%vx2l,& - & p%precv(level + 1)%map,info,work=work) + & p%precv(level + 1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -925,7 +937,8 @@ contains ! call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,& & zone,vy2l,& - & p%precv(level+1)%map,info,work=work) + & p%precv(level+1)%map,info,work=work,& + & vtx=wv(1),vty=wv(2)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -955,13 +968,13 @@ contains if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& & vty,zone,vy2l,& & base_desc, trans,& - & sweeps,work,info,init='Z') + & sweeps,work,wv,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') + & sweeps,work,wv,info,init='Z') end if if (info /= psb_success_) then @@ -1014,7 +1027,8 @@ contains 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) + & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& + & wv => p%precv(level)%wrk%wv) !Assemble rhs, w, v, v1, x diff --git a/mlprec/impl/mld_zprecaply.f90 b/mlprec/impl/mld_zprecaply.f90 index 80328261..2bb8a55c 100644 --- a/mlprec/impl/mld_zprecaply.f90 +++ b/mlprec/impl/mld_zprecaply.f90 @@ -323,6 +323,7 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work) complex(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_zprecaply' @@ -358,6 +359,10 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + if (size(prec%precv) >1) then ! ! Number of levels > 1: apply the multilevel preconditioner @@ -375,31 +380,29 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then - ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. - ! - twoside: block - type(psb_z_vect_type) :: w1,w2 - call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.) + associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,& + & wv => prec%precv(1)%wrk%wv) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! call psb_geaxpby(zone,x,zzero,w1,desc_data,info) select case(trans_) case ('N') do k=1, nswps call prec%precv(1)%sm%apply(zone,w1,zzero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm2a%apply(zone,w2,zzero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do - + case('T','C') do k=1, nswps call prec%precv(1)%sm2a%apply(zone,w1,zzero,w2,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) call prec%precv(1)%sm%apply(zone,w2,zzero,w1,desc_data,trans_,& - & ione, work_,info) + & ione, work_,wv,info) end do case default info = psb_err_from_subroutine_ @@ -407,14 +410,12 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work) goto 9999 end select call psb_geaxpby(zone,w1,zzero,y,desc_data,info) - call psb_gefree(w1,desc_data,info) - call psb_gefree(w2,desc_data,info) - end block twoside - - else - call prec%precv(1)%sm%apply(zone,x,zzero,y,desc_data,trans_,& - & nswps,work_,info) - end if + else + call prec%precv(1)%sm%apply(zone,x,zzero,y,desc_data,trans_,& + & nswps,work_,wv,info) + end if + end associate + else info = psb_err_from_subroutine_ai_ @@ -426,6 +427,7 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work) ! If the original distribution has an overlap we should fix that. call psb_halo(y,desc_data,info,data=psb_comm_mov_) + if (do_alloc_wrk) call prec%free_wrk(info) if (present(work)) then else @@ -459,10 +461,10 @@ subroutine mld_zprecaply1_vect(prec,x,desc_data,info,trans,work) ! Local variables character :: trans_ - type(psb_z_vect_type) :: ww complex(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act,iwsz, k, nswps + logical :: do_alloc_wrk character(len=20) :: name name='mld_zprecaply' @@ -498,63 +500,68 @@ subroutine mld_zprecaply1_vect(prec,x,desc_data,info,trans,work) call psb_errpush(info,name) goto 9999 end if - call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) - if (size(prec%precv) >1) then - ! - ! Number of levels > 1: apply the multilevel preconditioner - ! - call mld_mlprec_aply(zone,prec,x,zzero,ww,desc_data,trans_,work_,info) - if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_zmlprec_aply') - goto 9999 - end if + + do_alloc_wrk = .not.allocated(prec%wrk) + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) - else if (size(prec%precv) == 1) then - ! - ! Number of levels = 1: apply the base preconditioner - ! - nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) - if (allocated(prec%precv(1)%sm2a)) then + associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv) + call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.) + if (size(prec%precv) >1) then ! - ! This is a kludge for handling the symmetrized GS case. - ! Will need some rethinking. + ! Number of levels > 1: apply the multilevel preconditioner ! - select case(trans_) - case ('N') - do k=1, nswps - call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm2a%apply(zone,ww,zzero,x,desc_data,trans_,& - & ione, work_,info) - end do - case('T','C') - do k=1, nswps - call prec%precv(1)%sm2a%apply(zone,x,zzero,ww,desc_data,trans_,& - & ione, work_,info) - call prec%precv(1)%sm%apply(zone,ww,zzero,x,desc_data,trans_,& - & ione, work_,info) - end do - case default - info = psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='Invalid trans') - goto 9999 - end select - - else - call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,& - & nswps, work_,info) + call mld_mlprec_aply(zone,prec,x,zzero,ww,desc_data,trans_,work_,info) if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info) - end if - else + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_zmlprec_aply') + goto 9999 + end if - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='Invalid size of precv',& - & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) - goto 9999 - endif + else if (size(prec%precv) == 1) then + ! + ! Number of levels = 1: apply the base preconditioner + ! + nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) + if (allocated(prec%precv(1)%sm2a)) then + ! + ! This is a kludge for handling the symmetrized GS case. + ! Will need some rethinking. + ! + select case(trans_) + case ('N') + do k=1, nswps + call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm2a%apply(zone,ww,zzero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case('T','C') + do k=1, nswps + call prec%precv(1)%sm2a%apply(zone,x,zzero,ww,desc_data,trans_,& + & ione, work_,wv,info) + call prec%precv(1)%sm%apply(zone,ww,zzero,x,desc_data,trans_,& + & ione, work_,wv,info) + end do + case default + info = psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='Invalid trans') + goto 9999 + end select - if (info == 0) call psb_gefree(ww,desc_data,info) + else + call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,& + & nswps, work_,wv,info) + if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info) + end if + else + + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='Invalid size of precv',& + & i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) + goto 9999 + endif + end associate +!!$ if (info == 0) call psb_gefree(ww,desc_data,info) ! If the original distribution has an overlap we should fix that. call psb_halo(x,desc_data,info,data=psb_comm_mov_) diff --git a/mlprec/impl/mld_zprecset.F90 b/mlprec/impl/mld_zprecset.F90 index 4226b812..f5df9866 100644 --- a/mlprec/impl/mld_zprecset.F90 +++ b/mlprec/impl/mld_zprecset.F90 @@ -465,7 +465,7 @@ subroutine mld_zprecsetsm(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsm' info = psb_success_ @@ -528,7 +528,7 @@ subroutine mld_zprecsetsv(p,val,info,ilev,ilmax,pos) ! Local variables integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ - character(len=*), parameter :: name='mld_precseti' + character(len=*), parameter :: name='mld_precsetsv' info = psb_success_ diff --git a/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 index fa3d666a..6fc4ae38 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_c_vect_type),intent(inout), optional :: initu - type(psb_c_vect_type),intent(inout), optional :: wv(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i complex(psb_spk_), pointer :: aux(:) @@ -125,88 +125,95 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - ! Need to zero tx because of the apply_restr call. - call tx%zero() - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - if (info == 0) call sm%apply_restr(tx,trans_,aux,info) - if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) - - select case (init_) - case('Z') - call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(cone,y,czero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(cone,initu,czero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') + & a_err='invalid wv size in smoother_apply') goto 9999 - endif - - do i=1, sweeps-1 + end if + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + +!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) - if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit + + select case (init_) + case('Z') + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(cone,y,czero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do + + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + else info = psb_err_iarg_neg_ @@ -220,9 +227,9 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif - if (info ==0) call ww%free(info) - if (info ==0) call tx%free(info) - if (info ==0) call ty%free(info) +!!$ if (info ==0) call ww%free(info) +!!$ if (info ==0) call tx%free(info) +!!$ if (info ==0) call ty%free(info) if (info /= 0) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) diff --git a/mlprec/impl/smoother/mld_c_base_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_base_smoother_apply_vect.f90 index cc7eb484..9015b07c 100644 --- a/mlprec/impl/smoother/mld_c_base_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_base_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) use psb_base_mod use mld_c_base_smoother_mod, mld_protect_name => mld_c_base_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_c_vect_type),intent(inout), optional :: initu - type(psb_c_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: err_act character(len=20) :: name='c_base_smoother_apply' diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 index ce9d45f7..5e393745 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply_vect @@ -49,10 +49,10 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_c_vect_type),intent(inout), optional :: initu - type(psb_c_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: n_row,n_col type(psb_c_vect_type) :: tx, ty @@ -122,77 +122,86 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + if (size(wv) < 2) then + info = psb_err_internal_error_ + write(0,*) 'Size (WV) : ',size(wv) + call psb_errpush(info,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + associate(tx => wv(1), ty => wv(2)) +!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,y,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,initu,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - do i=1, sweeps-1 ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one AXPBY and one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + select case (init_) + case('Z') + + call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_geaxpby(cone,y,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - if (info /= psb_success_) exit + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit - if (info /= psb_success_) exit - end do + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info /= psb_success_) exit + end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - call tx%free(info) - if (info == psb_success_) call ty%free(info) - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='final cleanup with Jacobi sweeps > 1') - goto 9999 - end if + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if +!!$ call tx%free(info) +!!$ if (info == psb_success_) call ty%free(info) +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='final cleanup with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if + end associate + else info = psb_err_iarg_neg_ diff --git a/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 index 08784800..68bae713 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_d_vect_type),intent(inout), optional :: initu - type(psb_d_vect_type),intent(inout), optional :: wv(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i real(psb_dpk_), pointer :: aux(:) @@ -125,88 +125,95 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - ! Need to zero tx because of the apply_restr call. - call tx%zero() - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - if (info == 0) call sm%apply_restr(tx,trans_,aux,info) - if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) - - select case (init_) - case('Z') - call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(done,initu,dzero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') + & a_err='invalid wv size in smoother_apply') goto 9999 - endif - - do i=1, sweeps-1 + end if + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + +!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) - if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit + + select case (init_) + case('Z') + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(done,y,dzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do + + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + else info = psb_err_iarg_neg_ @@ -220,9 +227,9 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif - if (info ==0) call ww%free(info) - if (info ==0) call tx%free(info) - if (info ==0) call ty%free(info) +!!$ if (info ==0) call ww%free(info) +!!$ if (info ==0) call tx%free(info) +!!$ if (info ==0) call ty%free(info) if (info /= 0) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) diff --git a/mlprec/impl/smoother/mld_d_base_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_d_base_smoother_apply_vect.f90 index 2033b7a0..f153a18a 100644 --- a/mlprec/impl/smoother/mld_d_base_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_d_base_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) use psb_base_mod use mld_d_base_smoother_mod, mld_protect_name => mld_d_base_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_d_vect_type),intent(inout), optional :: initu - type(psb_d_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: err_act character(len=20) :: name='d_base_smoother_apply' diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 index 280bc233..c5d2b002 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply_vect @@ -49,10 +49,10 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_d_vect_type),intent(inout), optional :: initu - type(psb_d_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: n_row,n_col type(psb_d_vect_type) :: tx, ty @@ -122,77 +122,86 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + if (size(wv) < 2) then + info = psb_err_internal_error_ + write(0,*) 'Size (WV) : ',size(wv) + call psb_errpush(info,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + associate(tx => wv(1), ty => wv(2)) +!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,initu,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - do i=1, sweeps-1 ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one AXPBY and one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) + select case (init_) + case('Z') + + call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_geaxpby(done,y,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - if (info /= psb_success_) exit + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit - if (info /= psb_success_) exit - end do + call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info /= psb_success_) exit + end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - call tx%free(info) - if (info == psb_success_) call ty%free(info) - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='final cleanup with Jacobi sweeps > 1') - goto 9999 - end if + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if +!!$ call tx%free(info) +!!$ if (info == psb_success_) call ty%free(info) +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='final cleanup with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if + end associate + else info = psb_err_iarg_neg_ diff --git a/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 index 1e12bff5..ecfcd54b 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_s_vect_type),intent(inout), optional :: initu - type(psb_s_vect_type),intent(inout), optional :: wv(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i real(psb_spk_), pointer :: aux(:) @@ -125,88 +125,95 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - ! Need to zero tx because of the apply_restr call. - call tx%zero() - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - if (info == 0) call sm%apply_restr(tx,trans_,aux,info) - if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) - - select case (init_) - case('Z') - call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(sone,y,szero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(sone,initu,szero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') + & a_err='invalid wv size in smoother_apply') goto 9999 - endif - - do i=1, sweeps-1 + end if + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + +!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) - if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit + + select case (init_) + case('Z') + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(sone,y,szero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do + + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + else info = psb_err_iarg_neg_ @@ -220,9 +227,9 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif - if (info ==0) call ww%free(info) - if (info ==0) call tx%free(info) - if (info ==0) call ty%free(info) +!!$ if (info ==0) call ww%free(info) +!!$ if (info ==0) call tx%free(info) +!!$ if (info ==0) call ty%free(info) if (info /= 0) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) diff --git a/mlprec/impl/smoother/mld_s_base_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_base_smoother_apply_vect.f90 index 031b52e4..c3fb485e 100644 --- a/mlprec/impl/smoother/mld_s_base_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_base_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) use psb_base_mod use mld_s_base_smoother_mod, mld_protect_name => mld_s_base_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_s_vect_type),intent(inout), optional :: initu - type(psb_s_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: err_act character(len=20) :: name='s_base_smoother_apply' diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 index 24f81fa7..64927b4e 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply_vect @@ -49,10 +49,10 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_s_vect_type),intent(inout), optional :: initu - type(psb_s_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: n_row,n_col type(psb_s_vect_type) :: tx, ty @@ -122,77 +122,86 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + if (size(wv) < 2) then + info = psb_err_internal_error_ + write(0,*) 'Size (WV) : ',size(wv) + call psb_errpush(info,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + associate(tx => wv(1), ty => wv(2)) +!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,y,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,initu,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - do i=1, sweeps-1 ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one AXPBY and one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + select case (init_) + case('Z') + + call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_geaxpby(sone,y,szero,ty,desc_data,info) + call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - if (info /= psb_success_) exit + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit - if (info /= psb_success_) exit - end do + call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info /= psb_success_) exit + end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - call tx%free(info) - if (info == psb_success_) call ty%free(info) - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='final cleanup with Jacobi sweeps > 1') - goto 9999 - end if + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if +!!$ call tx%free(info) +!!$ if (info == psb_success_) call ty%free(info) +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='final cleanup with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if + end associate + else info = psb_err_iarg_neg_ diff --git a/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 index c0658c33..ae52866b 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_z_vect_type),intent(inout), optional :: initu - type(psb_z_vect_type),intent(inout), optional :: wv(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i complex(psb_dpk_), pointer :: aux(:) @@ -125,88 +125,95 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - ! Need to zero tx because of the apply_restr call. - call tx%zero() - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - if (info == 0) call sm%apply_restr(tx,trans_,aux,info) - if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) - - select case (init_) - case('Z') - call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) - if (info == 0) call sm%apply_restr(ty,trans_,aux,info) - if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') + & a_err='invalid wv size in smoother_apply') goto 9999 - endif - - do i=1, sweeps-1 + end if + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + +!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) - if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit + + select case (init_) + case('Z') + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select if (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do + + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + else info = psb_err_iarg_neg_ @@ -220,9 +227,9 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif - if (info ==0) call ww%free(info) - if (info ==0) call tx%free(info) - if (info ==0) call ty%free(info) +!!$ if (info ==0) call ww%free(info) +!!$ if (info ==0) call tx%free(info) +!!$ if (info ==0) call ty%free(info) if (info /= 0) then info = psb_err_alloc_dealloc_ call psb_errpush(info,name) diff --git a/mlprec/impl/smoother/mld_z_base_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_base_smoother_apply_vect.f90 index c8edb655..18e24b82 100644 --- a/mlprec/impl/smoother/mld_z_base_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_base_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) use psb_base_mod use mld_z_base_smoother_mod, mld_protect_name => mld_z_base_smoother_apply_vect implicit none @@ -48,10 +48,10 @@ subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_z_vect_type),intent(inout), optional :: initu - type(psb_z_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: err_act character(len=20) :: name='z_base_smoother_apply' diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 index 0e61fbab..af0e6488 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) use psb_base_mod use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply_vect @@ -49,10 +49,10 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_z_vect_type),intent(inout), optional :: initu - type(psb_z_vect_type),intent(inout), optional :: wv(:) ! integer(psb_ipk_) :: n_row,n_col type(psb_z_vect_type) :: tx, ty @@ -122,77 +122,86 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! to compute an approximate solution of a linear system. ! ! - call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) + if (size(wv) < 2) then + info = psb_err_internal_error_ + write(0,*) 'Size (WV) : ',size(wv) + call psb_errpush(info,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + associate(tx => wv(1), ty => wv(2)) +!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) +!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') - - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - do i=1, sweeps-1 ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one AXPBY and one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + select case (init_) + case('Z') + + call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') + + case('U') + if (.not.present(initu)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='missing initu to smoother_apply') + goto 9999 + end if + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - if (info /= psb_success_) exit + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit - if (info /= psb_success_) exit - end do + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info /= psb_success_) exit + end do - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - call tx%free(info) - if (info == psb_success_) call ty%free(info) - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='final cleanup with Jacobi sweeps > 1') - goto 9999 - end if + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if +!!$ call tx%free(info) +!!$ if (info == psb_success_) call ty%free(info) +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='final cleanup with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if + end associate + else info = psb_err_iarg_neg_ diff --git a/mlprec/mld_c_as_smoother.f90 b/mlprec/mld_c_as_smoother.f90 index c7d5ccd6..5e3f8f4e 100644 --- a/mlprec/mld_c_as_smoother.f90 +++ b/mlprec/mld_c_as_smoother.f90 @@ -181,7 +181,7 @@ module mld_c_as_smoother interface subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, & & psb_desc_type, psb_ipk_ @@ -194,10 +194,10 @@ module mld_c_as_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_c_vect_type),intent(inout), optional :: initu - type(psb_c_vect_type),intent(inout), optional :: wv(:) end subroutine mld_c_as_smoother_apply_vect end interface diff --git a/mlprec/mld_c_base_smoother_mod.f90 b/mlprec/mld_c_base_smoother_mod.f90 index 9c542cbd..b416bc79 100644 --- a/mlprec/mld_c_base_smoother_mod.f90 +++ b/mlprec/mld_c_base_smoother_mod.f90 @@ -158,7 +158,7 @@ module mld_c_base_smoother_mod interface subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, & & mld_c_base_smoother_type, psb_ipk_ @@ -170,10 +170,10 @@ module mld_c_base_smoother_mod character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_c_vect_type),intent(inout), optional :: initu - type(psb_c_vect_type),intent(inout), optional :: wv(:) end subroutine mld_c_base_smoother_apply_vect end interface diff --git a/mlprec/mld_c_jac_smoother.f90 b/mlprec/mld_c_jac_smoother.f90 index 3b464b6d..1c2221bb 100644 --- a/mlprec/mld_c_jac_smoother.f90 +++ b/mlprec/mld_c_jac_smoother.f90 @@ -85,7 +85,7 @@ module mld_c_jac_smoother interface subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, & & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,& & psb_ipk_ @@ -98,10 +98,10 @@ module mld_c_jac_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_c_vect_type),intent(inout), optional :: initu - type(psb_c_vect_type),intent(inout), optional :: wv(:) end subroutine mld_c_jac_smoother_apply_vect end interface diff --git a/mlprec/mld_c_onelev_mod.f90 b/mlprec/mld_c_onelev_mod.f90 index b76119a3..b4a69b53 100644 --- a/mlprec/mld_c_onelev_mod.f90 +++ b/mlprec/mld_c_onelev_mod.f90 @@ -583,10 +583,9 @@ contains integer(psb_ipk_), intent(out) :: info class(psb_c_base_vect_type), intent(in), optional :: vmold ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv, i 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) @@ -599,7 +598,12 @@ contains call psb_geasb(lv%wrk%vty,& & lv%base_desc,info,& & scratch=.true.,mold=vmold) - + allocate(lv%wrk%wv(nwv),stat=info) + do i=1,nwv + call psb_geasb(lv%wrk%wv(i),& + & lv%base_desc,info,& + & scratch=.true.,mold=vmold) + end do end subroutine c_base_onelev_allocate_wrk @@ -609,12 +613,17 @@ contains class(mld_c_onelev_type), target, intent(inout) :: lv integer(psb_ipk_), intent(out) :: info ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv,i 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) + nwv = size(lv%wrk%wv) + do i=1,nwv + call lv%wrk%wv(i)%free(info) + end do + end subroutine c_base_onelev_free_wrk end module mld_c_onelev_mod diff --git a/mlprec/mld_d_as_smoother.f90 b/mlprec/mld_d_as_smoother.f90 index 18d2313c..4ffbbde6 100644 --- a/mlprec/mld_d_as_smoother.f90 +++ b/mlprec/mld_d_as_smoother.f90 @@ -181,7 +181,7 @@ module mld_d_as_smoother interface subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, & & psb_desc_type, psb_ipk_ @@ -194,10 +194,10 @@ module mld_d_as_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_d_vect_type),intent(inout), optional :: initu - type(psb_d_vect_type),intent(inout), optional :: wv(:) end subroutine mld_d_as_smoother_apply_vect end interface diff --git a/mlprec/mld_d_base_smoother_mod.f90 b/mlprec/mld_d_base_smoother_mod.f90 index b42c96e6..e62fbc16 100644 --- a/mlprec/mld_d_base_smoother_mod.f90 +++ b/mlprec/mld_d_base_smoother_mod.f90 @@ -158,7 +158,7 @@ module mld_d_base_smoother_mod interface subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, & & mld_d_base_smoother_type, psb_ipk_ @@ -170,10 +170,10 @@ module mld_d_base_smoother_mod character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_d_vect_type),intent(inout), optional :: initu - type(psb_d_vect_type),intent(inout), optional :: wv(:) end subroutine mld_d_base_smoother_apply_vect end interface diff --git a/mlprec/mld_d_jac_smoother.f90 b/mlprec/mld_d_jac_smoother.f90 index ef03d431..be03b750 100644 --- a/mlprec/mld_d_jac_smoother.f90 +++ b/mlprec/mld_d_jac_smoother.f90 @@ -85,7 +85,7 @@ module mld_d_jac_smoother interface subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) import :: psb_desc_type, mld_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, & & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,& & psb_ipk_ @@ -98,10 +98,10 @@ module mld_d_jac_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) + type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_d_vect_type),intent(inout), optional :: initu - type(psb_d_vect_type),intent(inout), optional :: wv(:) end subroutine mld_d_jac_smoother_apply_vect end interface diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index ce9d0157..636d566c 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -583,10 +583,9 @@ contains integer(psb_ipk_), intent(out) :: info class(psb_d_base_vect_type), intent(in), optional :: vmold ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv, i 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) @@ -599,7 +598,12 @@ contains call psb_geasb(lv%wrk%vty,& & lv%base_desc,info,& & scratch=.true.,mold=vmold) - + allocate(lv%wrk%wv(nwv),stat=info) + do i=1,nwv + call psb_geasb(lv%wrk%wv(i),& + & lv%base_desc,info,& + & scratch=.true.,mold=vmold) + end do end subroutine d_base_onelev_allocate_wrk @@ -609,12 +613,17 @@ contains class(mld_d_onelev_type), target, intent(inout) :: lv integer(psb_ipk_), intent(out) :: info ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv,i 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) + nwv = size(lv%wrk%wv) + do i=1,nwv + call lv%wrk%wv(i)%free(info) + end do + end subroutine d_base_onelev_free_wrk end module mld_d_onelev_mod diff --git a/mlprec/mld_s_as_smoother.f90 b/mlprec/mld_s_as_smoother.f90 index e800c279..eab3f721 100644 --- a/mlprec/mld_s_as_smoother.f90 +++ b/mlprec/mld_s_as_smoother.f90 @@ -181,7 +181,7 @@ module mld_s_as_smoother interface subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & & psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, & & psb_desc_type, psb_ipk_ @@ -194,10 +194,10 @@ module mld_s_as_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_s_vect_type),intent(inout), optional :: initu - type(psb_s_vect_type),intent(inout), optional :: wv(:) end subroutine mld_s_as_smoother_apply_vect end interface diff --git a/mlprec/mld_s_base_smoother_mod.f90 b/mlprec/mld_s_base_smoother_mod.f90 index b5cb6020..dcf3188a 100644 --- a/mlprec/mld_s_base_smoother_mod.f90 +++ b/mlprec/mld_s_base_smoother_mod.f90 @@ -158,7 +158,7 @@ module mld_s_base_smoother_mod interface subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, & & mld_s_base_smoother_type, psb_ipk_ @@ -170,10 +170,10 @@ module mld_s_base_smoother_mod character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_s_vect_type),intent(inout), optional :: initu - type(psb_s_vect_type),intent(inout), optional :: wv(:) end subroutine mld_s_base_smoother_apply_vect end interface diff --git a/mlprec/mld_s_jac_smoother.f90 b/mlprec/mld_s_jac_smoother.f90 index ccb7d896..b13f9cbe 100644 --- a/mlprec/mld_s_jac_smoother.f90 +++ b/mlprec/mld_s_jac_smoother.f90 @@ -85,7 +85,7 @@ module mld_s_jac_smoother interface subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) import :: psb_desc_type, mld_s_jac_smoother_type, psb_s_vect_type, psb_spk_, & & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,& & psb_ipk_ @@ -98,10 +98,10 @@ module mld_s_jac_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_s_vect_type),intent(inout), optional :: initu - type(psb_s_vect_type),intent(inout), optional :: wv(:) end subroutine mld_s_jac_smoother_apply_vect end interface diff --git a/mlprec/mld_s_onelev_mod.f90 b/mlprec/mld_s_onelev_mod.f90 index 88f6bf14..94071867 100644 --- a/mlprec/mld_s_onelev_mod.f90 +++ b/mlprec/mld_s_onelev_mod.f90 @@ -583,10 +583,9 @@ contains integer(psb_ipk_), intent(out) :: info class(psb_s_base_vect_type), intent(in), optional :: vmold ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv, i 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) @@ -599,7 +598,12 @@ contains call psb_geasb(lv%wrk%vty,& & lv%base_desc,info,& & scratch=.true.,mold=vmold) - + allocate(lv%wrk%wv(nwv),stat=info) + do i=1,nwv + call psb_geasb(lv%wrk%wv(i),& + & lv%base_desc,info,& + & scratch=.true.,mold=vmold) + end do end subroutine s_base_onelev_allocate_wrk @@ -609,12 +613,17 @@ contains class(mld_s_onelev_type), target, intent(inout) :: lv integer(psb_ipk_), intent(out) :: info ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv,i 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) + nwv = size(lv%wrk%wv) + do i=1,nwv + call lv%wrk%wv(i)%free(info) + end do + end subroutine s_base_onelev_free_wrk end module mld_s_onelev_mod diff --git a/mlprec/mld_z_as_smoother.f90 b/mlprec/mld_z_as_smoother.f90 index 468026c3..9254ba17 100644 --- a/mlprec/mld_z_as_smoother.f90 +++ b/mlprec/mld_z_as_smoother.f90 @@ -181,7 +181,7 @@ module mld_z_as_smoother interface subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & & psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, & & psb_desc_type, psb_ipk_ @@ -194,10 +194,10 @@ module mld_z_as_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_z_vect_type),intent(inout), optional :: initu - type(psb_z_vect_type),intent(inout), optional :: wv(:) end subroutine mld_z_as_smoother_apply_vect end interface diff --git a/mlprec/mld_z_base_smoother_mod.f90 b/mlprec/mld_z_base_smoother_mod.f90 index 9b607ecf..505b4926 100644 --- a/mlprec/mld_z_base_smoother_mod.f90 +++ b/mlprec/mld_z_base_smoother_mod.f90 @@ -158,7 +158,7 @@ module mld_z_base_smoother_mod interface subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& - & trans,sweeps,work,info,init,initu,wv) + & trans,sweeps,work,wv,info,init,initu) import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, & & mld_z_base_smoother_type, psb_ipk_ @@ -170,10 +170,10 @@ module mld_z_base_smoother_mod character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_z_vect_type),intent(inout), optional :: initu - type(psb_z_vect_type),intent(inout), optional :: wv(:) end subroutine mld_z_base_smoother_apply_vect end interface diff --git a/mlprec/mld_z_jac_smoother.f90 b/mlprec/mld_z_jac_smoother.f90 index 08cf38f3..07b30681 100644 --- a/mlprec/mld_z_jac_smoother.f90 +++ b/mlprec/mld_z_jac_smoother.f90 @@ -85,7 +85,7 @@ module mld_z_jac_smoother interface subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& - & sweeps,work,info,init,initu,wv) + & sweeps,work,wv,info,init,initu) import :: psb_desc_type, mld_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, & & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,& & psb_ipk_ @@ -98,10 +98,10 @@ module mld_z_jac_smoother character(len=1),intent(in) :: trans integer(psb_ipk_), intent(in) :: sweeps complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info character, intent(in), optional :: init type(psb_z_vect_type),intent(inout), optional :: initu - type(psb_z_vect_type),intent(inout), optional :: wv(:) end subroutine mld_z_jac_smoother_apply_vect end interface diff --git a/mlprec/mld_z_onelev_mod.f90 b/mlprec/mld_z_onelev_mod.f90 index 94a8fe43..4476e616 100644 --- a/mlprec/mld_z_onelev_mod.f90 +++ b/mlprec/mld_z_onelev_mod.f90 @@ -583,10 +583,9 @@ contains integer(psb_ipk_), intent(out) :: info class(psb_z_base_vect_type), intent(in), optional :: vmold ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv, i 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) @@ -599,7 +598,12 @@ contains call psb_geasb(lv%wrk%vty,& & lv%base_desc,info,& & scratch=.true.,mold=vmold) - + allocate(lv%wrk%wv(nwv),stat=info) + do i=1,nwv + call psb_geasb(lv%wrk%wv(i),& + & lv%base_desc,info,& + & scratch=.true.,mold=vmold) + end do end subroutine z_base_onelev_allocate_wrk @@ -609,12 +613,17 @@ contains class(mld_z_onelev_type), target, intent(inout) :: lv integer(psb_ipk_), intent(out) :: info ! - integer(psb_ipk_) :: nwv + integer(psb_ipk_) :: nwv,i 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) + nwv = size(lv%wrk%wv) + do i=1,nwv + call lv%wrk%wv(i)%free(info) + end do + end subroutine z_base_onelev_free_wrk end module mld_z_onelev_mod