From 28620a78798a691effecdfadd5f74a9e2a68e593 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 3 May 2019 16:12:12 +0100 Subject: [PATCH] New wrk methods in psblas%prec --- krylov/psb_ckrylov.f90 | 7 +- krylov/psb_dkrylov.f90 | 7 +- krylov/psb_skrylov.f90 | 7 +- krylov/psb_zkrylov.f90 | 7 +- prec/impl/psb_c_bjacprec_impl.f90 | 93 +++++++++++++------------ prec/impl/psb_d_bjacprec_impl.f90 | 93 +++++++++++++------------ prec/impl/psb_s_bjacprec_impl.f90 | 93 +++++++++++++------------ prec/impl/psb_z_bjacprec_impl.f90 | 93 +++++++++++++------------ prec/psb_c_base_prec_mod.f90 | 84 ++++++++++++++++++++++ prec/psb_c_bjacprec.f90 | 112 +++++++++++++++++++++++++++++- prec/psb_c_prec_type.f90 | 87 +++++++++++++++++++++++ prec/psb_d_base_prec_mod.f90 | 84 ++++++++++++++++++++++ prec/psb_d_bjacprec.f90 | 112 +++++++++++++++++++++++++++++- prec/psb_d_prec_type.f90 | 87 +++++++++++++++++++++++ prec/psb_s_base_prec_mod.f90 | 84 ++++++++++++++++++++++ prec/psb_s_bjacprec.f90 | 112 +++++++++++++++++++++++++++++- prec/psb_s_prec_type.f90 | 87 +++++++++++++++++++++++ prec/psb_z_base_prec_mod.f90 | 84 ++++++++++++++++++++++ prec/psb_z_bjacprec.f90 | 112 +++++++++++++++++++++++++++++- prec/psb_z_prec_type.f90 | 87 +++++++++++++++++++++++ 20 files changed, 1344 insertions(+), 188 deletions(-) diff --git a/krylov/psb_ckrylov.f90 b/krylov/psb_ckrylov.f90 index e5c36078..2aae7913 100644 --- a/krylov/psb_ckrylov.f90 +++ b/krylov/psb_ckrylov.f90 @@ -151,7 +151,7 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& procedure(psb_ckryl_rest_vect) :: psb_crgmres_vect, psb_ccgstabl_vect, psb_cgcr_vect procedure(psb_ckryl_cond_vect) :: psb_ccg_vect, psb_cfcg_vect - + logical :: do_alloc_wrk integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_ character(len=20) :: name @@ -172,6 +172,9 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& itrace_ = -1 end if + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + select case(psb_toupper(method)) case('CG') call psb_ccg_vect(a,prec,b,x,eps,desc_a,info,& @@ -205,6 +208,8 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,& &itmax,iter,err,itrace=itrace_,istop=istop) end select + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + if(info /= psb_success_) then call psb_errpush(info,name) goto 9999 diff --git a/krylov/psb_dkrylov.f90 b/krylov/psb_dkrylov.f90 index c4049f73..1f887299 100644 --- a/krylov/psb_dkrylov.f90 +++ b/krylov/psb_dkrylov.f90 @@ -151,7 +151,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& procedure(psb_dkryl_rest_vect) :: psb_drgmres_vect, psb_dcgstabl_vect, psb_dgcr_vect procedure(psb_dkryl_cond_vect) :: psb_dcg_vect, psb_dfcg_vect - + logical :: do_alloc_wrk integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_ character(len=20) :: name @@ -172,6 +172,9 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& itrace_ = -1 end if + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + select case(psb_toupper(method)) case('CG') call psb_dcg_vect(a,prec,b,x,eps,desc_a,info,& @@ -205,6 +208,8 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& &itmax,iter,err,itrace=itrace_,istop=istop) end select + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + if(info /= psb_success_) then call psb_errpush(info,name) goto 9999 diff --git a/krylov/psb_skrylov.f90 b/krylov/psb_skrylov.f90 index 72a6eae8..2dc20812 100644 --- a/krylov/psb_skrylov.f90 +++ b/krylov/psb_skrylov.f90 @@ -151,7 +151,7 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& procedure(psb_skryl_rest_vect) :: psb_srgmres_vect, psb_scgstabl_vect, psb_sgcr_vect procedure(psb_skryl_cond_vect) :: psb_scg_vect, psb_sfcg_vect - + logical :: do_alloc_wrk integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_ character(len=20) :: name @@ -172,6 +172,9 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& itrace_ = -1 end if + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + select case(psb_toupper(method)) case('CG') call psb_scg_vect(a,prec,b,x,eps,desc_a,info,& @@ -205,6 +208,8 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,& &itmax,iter,err,itrace=itrace_,istop=istop) end select + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + if(info /= psb_success_) then call psb_errpush(info,name) goto 9999 diff --git a/krylov/psb_zkrylov.f90 b/krylov/psb_zkrylov.f90 index 2fbb5ba7..626499df 100644 --- a/krylov/psb_zkrylov.f90 +++ b/krylov/psb_zkrylov.f90 @@ -151,7 +151,7 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& procedure(psb_zkryl_rest_vect) :: psb_zrgmres_vect, psb_zcgstabl_vect, psb_zgcr_vect procedure(psb_zkryl_cond_vect) :: psb_zcg_vect, psb_zfcg_vect - + logical :: do_alloc_wrk integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_ character(len=20) :: name @@ -172,6 +172,9 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& itrace_ = -1 end if + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a) + select case(psb_toupper(method)) case('CG') call psb_zcg_vect(a,prec,b,x,eps,desc_a,info,& @@ -205,6 +208,8 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,& &itmax,iter,err,itrace=itrace_,istop=istop) end select + if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info) + if(info /= psb_success_) then call psb_errpush(info,name) goto 9999 diff --git a/prec/impl/psb_c_bjacprec_impl.f90 b/prec/impl/psb_c_bjacprec_impl.f90 index c1e24f42..de453684 100644 --- a/prec/impl/psb_c_bjacprec_impl.f90 +++ b/prec/impl/psb_c_bjacprec_impl.f90 @@ -90,6 +90,7 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit + logical :: do_alloc_wrk character :: trans_ character(len=20) :: name='c_bjac_prec_apply' character(len=20) :: ch_err @@ -154,55 +155,57 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) - - select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - - select case(trans_) - case('N') - call psb_spsm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - case('T') - call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - case('C') - - call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - call wv1%mlt(cone,prec%dv,wv,czero,info,conjgx=trans_) - - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - end select - if (info /= psb_success_) then - ch_err="psb_spsm" + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + + associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) + + select case(prec%iprcparm(psb_f_type_)) + case(psb_f_ilu_n_) + + select case(trans_) + case('N') + call psb_spsm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + case('T') + call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case('C') + + call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + call wv1%mlt(cone,prec%dv,wv,czero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid factorization') goto 9999 - end if - - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Invalid factorization') - goto 9999 - end select + end select + end associate + call psb_halo(y,desc_data,info,data=psb_comm_mov_) - call wv%free(info) - call wv1%free(info) + if (do_alloc_wrk) call prec%free_wrk(info) if (n_col <= size(work)) then if ((4*n_col+n_col) <= size(work)) then else diff --git a/prec/impl/psb_d_bjacprec_impl.f90 b/prec/impl/psb_d_bjacprec_impl.f90 index 8420da45..ef8c52c3 100644 --- a/prec/impl/psb_d_bjacprec_impl.f90 +++ b/prec/impl/psb_d_bjacprec_impl.f90 @@ -90,6 +90,7 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit + logical :: do_alloc_wrk character :: trans_ character(len=20) :: name='d_bjac_prec_apply' character(len=20) :: ch_err @@ -154,55 +155,57 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) - - select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - - select case(trans_) - case('N') - call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - case('T') - call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - case('C') - - call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - call wv1%mlt(done,prec%dv,wv,dzero,info,conjgx=trans_) - - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - end select - if (info /= psb_success_) then - ch_err="psb_spsm" + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + + associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) + + select case(prec%iprcparm(psb_f_type_)) + case(psb_f_ilu_n_) + + select case(trans_) + case('N') + call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + case('T') + call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case('C') + + call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + call wv1%mlt(done,prec%dv,wv,dzero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid factorization') goto 9999 - end if - - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Invalid factorization') - goto 9999 - end select + end select + end associate + call psb_halo(y,desc_data,info,data=psb_comm_mov_) - call wv%free(info) - call wv1%free(info) + if (do_alloc_wrk) call prec%free_wrk(info) if (n_col <= size(work)) then if ((4*n_col+n_col) <= size(work)) then else diff --git a/prec/impl/psb_s_bjacprec_impl.f90 b/prec/impl/psb_s_bjacprec_impl.f90 index eecf85d1..3a9cfce2 100644 --- a/prec/impl/psb_s_bjacprec_impl.f90 +++ b/prec/impl/psb_s_bjacprec_impl.f90 @@ -90,6 +90,7 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit + logical :: do_alloc_wrk character :: trans_ character(len=20) :: name='s_bjac_prec_apply' character(len=20) :: ch_err @@ -154,55 +155,57 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) - - select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - - select case(trans_) - case('N') - call psb_spsm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - case('T') - call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - case('C') - - call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - call wv1%mlt(sone,prec%dv,wv,szero,info,conjgx=trans_) - - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - end select - if (info /= psb_success_) then - ch_err="psb_spsm" + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + + associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) + + select case(prec%iprcparm(psb_f_type_)) + case(psb_f_ilu_n_) + + select case(trans_) + case('N') + call psb_spsm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + case('T') + call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case('C') + + call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + call wv1%mlt(sone,prec%dv,wv,szero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid factorization') goto 9999 - end if - - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Invalid factorization') - goto 9999 - end select + end select + end associate + call psb_halo(y,desc_data,info,data=psb_comm_mov_) - call wv%free(info) - call wv1%free(info) + if (do_alloc_wrk) call prec%free_wrk(info) if (n_col <= size(work)) then if ((4*n_col+n_col) <= size(work)) then else diff --git a/prec/impl/psb_z_bjacprec_impl.f90 b/prec/impl/psb_z_bjacprec_impl.f90 index 57825bf4..b70018f4 100644 --- a/prec/impl/psb_z_bjacprec_impl.f90 +++ b/prec/impl/psb_z_bjacprec_impl.f90 @@ -90,6 +90,7 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act, ierr(5) integer(psb_ipk_) :: debug_level, debug_unit + logical :: do_alloc_wrk character :: trans_ character(len=20) :: name='z_bjac_prec_apply' character(len=20) :: ch_err @@ -154,55 +155,57 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) - - select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - - select case(trans_) - case('N') - call psb_spsm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - case('T') - call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& - & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - case('C') - - call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_, work=aux) - - call wv1%mlt(zone,prec%dv,wv,zzero,info,conjgx=trans_) - - if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& - & beta,y,desc_data,info,& - & trans=trans_,scale='U',choice=psb_none_,work=aux) - - end select - if (info /= psb_success_) then - ch_err="psb_spsm" + do_alloc_wrk = .not.prec%is_allocated_wrk() + if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) + + associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) + + select case(prec%iprcparm(psb_f_type_)) + case(psb_f_ilu_n_) + + select case(trans_) + case('N') + call psb_spsm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + case('T') + call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case('C') + + call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + call wv1%mlt(zone,prec%dv,wv,zzero,info,conjgx=trans_) + + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& + & beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + end select + if (info /= psb_success_) then + ch_err="psb_spsm" + goto 9999 + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid factorization') goto 9999 - end if - - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='Invalid factorization') - goto 9999 - end select + end select + end associate + call psb_halo(y,desc_data,info,data=psb_comm_mov_) - call wv%free(info) - call wv1%free(info) + if (do_alloc_wrk) call prec%free_wrk(info) if (n_col <= size(work)) then if ((4*n_col+n_col) <= size(work)) then else diff --git a/prec/psb_c_base_prec_mod.f90 b/prec/psb_c_base_prec_mod.f90 index 9fe9d8e1..8b453241 100644 --- a/prec/psb_c_base_prec_mod.f90 +++ b/prec/psb_c_base_prec_mod.f90 @@ -61,6 +61,9 @@ module psb_c_base_prec_mod generic, public :: build => precbld generic, public :: descr => precdescr procedure, pass(prec) :: desc_prefix => psb_c_base_desc_prefix + procedure, pass(prec) :: allocate_wrk => psb_c_base_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_c_base_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_c_base_is_allocated_wrk procedure(psb_c_base_precbld), pass(prec), deferred :: precbld procedure(psb_c_base_sizeof), pass(prec), deferred :: sizeof procedure(psb_c_base_precinit), pass(prec), deferred :: precinit @@ -258,6 +261,87 @@ contains end subroutine psb_c_base_precsetc + subroutine psb_c_base_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_c_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_c_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_c_base_allocate_wrk + + subroutine psb_c_base_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_c_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_c_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_c_base_free_wrk + + function psb_c_base_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_c_base_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = .true. + + end function psb_c_base_is_allocated_wrk + subroutine psb_c_base_set_ctxt(prec,ictxt) implicit none class(psb_c_base_prec_type), intent(inout) :: prec diff --git a/prec/psb_c_bjacprec.f90 b/prec/psb_c_bjacprec.f90 index 7da40e38..58676962 100644 --- a/prec/psb_c_bjacprec.f90 +++ b/prec/psb_c_bjacprec.f90 @@ -36,7 +36,7 @@ module psb_c_bjacprec type, extends(psb_c_base_prec_type) :: psb_c_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) type(psb_cspmat_type), allocatable :: av(:) - type(psb_c_vect_type), allocatable :: dv + type(psb_c_vect_type), allocatable :: dv, wrk(:) contains procedure, pass(prec) :: c_apply_v => psb_c_bjac_apply_vect procedure, pass(prec) :: c_apply => psb_c_bjac_apply @@ -49,6 +49,9 @@ module psb_c_bjacprec procedure, pass(prec) :: free => psb_c_bjac_precfree procedure, pass(prec) :: sizeof => psb_c_bjac_sizeof procedure, pass(prec) :: get_nzeros => psb_c_bjac_get_nzeros + procedure, pass(prec) :: allocate_wrk => psb_c_bjac_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_c_bjac_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_c_bjac_is_allocated_wrk end type psb_c_bjac_prec_type private :: psb_c_bjac_sizeof, psb_c_bjac_precdescr, psb_c_bjac_get_nzeros @@ -308,4 +311,111 @@ contains end subroutine psb_c_bjac_clone + subroutine psb_c_bjac_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_c_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act, i + character(len=20) :: name + + info=psb_success_ + name = 'psb_c_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + if (allocated(prec%wrk)) then + if (size(prec%wrk)<2) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + if (.not.allocated(prec%wrk)) then + if (.not.present(desc)) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 + end if + allocate(prec%wrk(2),stat=info) + do i=1, 2 + if (info == 0) call psb_geall(prec%wrk(i),desc,info) + if (info == 0) call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.) + end do + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_c_bjac_allocate_wrk + + subroutine psb_c_bjac_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_c_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: i + character(len=20) :: name + + info=psb_success_ + name = 'psb_c_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + info = psb_success_ + if (allocated(prec%wrk)) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_c_bjac_free_wrk + + function psb_c_bjac_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_c_bjac_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = allocated(prec%wrk) + + end function psb_c_bjac_is_allocated_wrk + + end module psb_c_bjacprec diff --git a/prec/psb_c_prec_type.f90 b/prec/psb_c_prec_type.f90 index 5e341e1c..6d9aa908 100644 --- a/prec/psb_c_prec_type.f90 +++ b/prec/psb_c_prec_type.f90 @@ -54,6 +54,9 @@ module psb_c_prec_type procedure, pass(prec) :: build => psb_cprecbld procedure, pass(prec) :: init => psb_cprecinit procedure, pass(prec) :: descr => psb_cfile_prec_descr + procedure, pass(prec) :: allocate_wrk => psb_c_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_c_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_c_is_allocated_wrk end type psb_cprec_type interface psb_precfree @@ -193,6 +196,90 @@ contains end subroutine psb_c_prec_dump + subroutine psb_c_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_cprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_c_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner' + return + end if + + call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_c_allocate_wrk + + subroutine psb_c_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_cprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_c_free_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to free a non-built preconditioner' + return + end if + + call prec%prec%free_wrk(info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_c_free_wrk + + function psb_c_is_allocated_wrk(prec) result(res) + implicit none + + ! Arguments + class(psb_cprec_type), intent(in) :: prec + logical :: res + + if (.not.allocated(prec%prec)) then + res = .false. + else + res = prec%prec%is_allocated_wrk() + end if + + end function psb_c_is_allocated_wrk subroutine psb_c_precfree(p,info) use psb_base_mod diff --git a/prec/psb_d_base_prec_mod.f90 b/prec/psb_d_base_prec_mod.f90 index 738f6573..dc39ca70 100644 --- a/prec/psb_d_base_prec_mod.f90 +++ b/prec/psb_d_base_prec_mod.f90 @@ -61,6 +61,9 @@ module psb_d_base_prec_mod generic, public :: build => precbld generic, public :: descr => precdescr procedure, pass(prec) :: desc_prefix => psb_d_base_desc_prefix + procedure, pass(prec) :: allocate_wrk => psb_d_base_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_d_base_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_d_base_is_allocated_wrk procedure(psb_d_base_precbld), pass(prec), deferred :: precbld procedure(psb_d_base_sizeof), pass(prec), deferred :: sizeof procedure(psb_d_base_precinit), pass(prec), deferred :: precinit @@ -258,6 +261,87 @@ contains end subroutine psb_d_base_precsetc + subroutine psb_d_base_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_d_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_d_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_d_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_d_base_allocate_wrk + + subroutine psb_d_base_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_d_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_d_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_d_base_free_wrk + + function psb_d_base_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_d_base_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = .true. + + end function psb_d_base_is_allocated_wrk + subroutine psb_d_base_set_ctxt(prec,ictxt) implicit none class(psb_d_base_prec_type), intent(inout) :: prec diff --git a/prec/psb_d_bjacprec.f90 b/prec/psb_d_bjacprec.f90 index 3568af07..5e7413b8 100644 --- a/prec/psb_d_bjacprec.f90 +++ b/prec/psb_d_bjacprec.f90 @@ -36,7 +36,7 @@ module psb_d_bjacprec type, extends(psb_d_base_prec_type) :: psb_d_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) type(psb_dspmat_type), allocatable :: av(:) - type(psb_d_vect_type), allocatable :: dv + type(psb_d_vect_type), allocatable :: dv, wrk(:) contains procedure, pass(prec) :: d_apply_v => psb_d_bjac_apply_vect procedure, pass(prec) :: d_apply => psb_d_bjac_apply @@ -49,6 +49,9 @@ module psb_d_bjacprec procedure, pass(prec) :: free => psb_d_bjac_precfree procedure, pass(prec) :: sizeof => psb_d_bjac_sizeof procedure, pass(prec) :: get_nzeros => psb_d_bjac_get_nzeros + procedure, pass(prec) :: allocate_wrk => psb_d_bjac_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_d_bjac_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_d_bjac_is_allocated_wrk end type psb_d_bjac_prec_type private :: psb_d_bjac_sizeof, psb_d_bjac_precdescr, psb_d_bjac_get_nzeros @@ -308,4 +311,111 @@ contains end subroutine psb_d_bjac_clone + subroutine psb_d_bjac_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_d_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_d_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act, i + character(len=20) :: name + + info=psb_success_ + name = 'psb_d_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + if (allocated(prec%wrk)) then + if (size(prec%wrk)<2) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + if (.not.allocated(prec%wrk)) then + if (.not.present(desc)) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 + end if + allocate(prec%wrk(2),stat=info) + do i=1, 2 + if (info == 0) call psb_geall(prec%wrk(i),desc,info) + if (info == 0) call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.) + end do + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_d_bjac_allocate_wrk + + subroutine psb_d_bjac_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_d_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: i + character(len=20) :: name + + info=psb_success_ + name = 'psb_d_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + info = psb_success_ + if (allocated(prec%wrk)) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_d_bjac_free_wrk + + function psb_d_bjac_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_d_bjac_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = allocated(prec%wrk) + + end function psb_d_bjac_is_allocated_wrk + + end module psb_d_bjacprec diff --git a/prec/psb_d_prec_type.f90 b/prec/psb_d_prec_type.f90 index f420d282..f50339a3 100644 --- a/prec/psb_d_prec_type.f90 +++ b/prec/psb_d_prec_type.f90 @@ -54,6 +54,9 @@ module psb_d_prec_type procedure, pass(prec) :: build => psb_dprecbld procedure, pass(prec) :: init => psb_dprecinit procedure, pass(prec) :: descr => psb_dfile_prec_descr + procedure, pass(prec) :: allocate_wrk => psb_d_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_d_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_d_is_allocated_wrk end type psb_dprec_type interface psb_precfree @@ -193,6 +196,90 @@ contains end subroutine psb_d_prec_dump + subroutine psb_d_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_dprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_d_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_d_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner' + return + end if + + call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_d_allocate_wrk + + subroutine psb_d_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_dprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_d_free_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to free a non-built preconditioner' + return + end if + + call prec%prec%free_wrk(info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_d_free_wrk + + function psb_d_is_allocated_wrk(prec) result(res) + implicit none + + ! Arguments + class(psb_dprec_type), intent(in) :: prec + logical :: res + + if (.not.allocated(prec%prec)) then + res = .false. + else + res = prec%prec%is_allocated_wrk() + end if + + end function psb_d_is_allocated_wrk subroutine psb_d_precfree(p,info) use psb_base_mod diff --git a/prec/psb_s_base_prec_mod.f90 b/prec/psb_s_base_prec_mod.f90 index eb468301..b2a89f06 100644 --- a/prec/psb_s_base_prec_mod.f90 +++ b/prec/psb_s_base_prec_mod.f90 @@ -61,6 +61,9 @@ module psb_s_base_prec_mod generic, public :: build => precbld generic, public :: descr => precdescr procedure, pass(prec) :: desc_prefix => psb_s_base_desc_prefix + procedure, pass(prec) :: allocate_wrk => psb_s_base_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_s_base_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_s_base_is_allocated_wrk procedure(psb_s_base_precbld), pass(prec), deferred :: precbld procedure(psb_s_base_sizeof), pass(prec), deferred :: sizeof procedure(psb_s_base_precinit), pass(prec), deferred :: precinit @@ -258,6 +261,87 @@ contains end subroutine psb_s_base_precsetc + subroutine psb_s_base_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_s_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_s_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_s_base_allocate_wrk + + subroutine psb_s_base_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_s_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_s_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_s_base_free_wrk + + function psb_s_base_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_s_base_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = .true. + + end function psb_s_base_is_allocated_wrk + subroutine psb_s_base_set_ctxt(prec,ictxt) implicit none class(psb_s_base_prec_type), intent(inout) :: prec diff --git a/prec/psb_s_bjacprec.f90 b/prec/psb_s_bjacprec.f90 index 093da0ca..aba98871 100644 --- a/prec/psb_s_bjacprec.f90 +++ b/prec/psb_s_bjacprec.f90 @@ -36,7 +36,7 @@ module psb_s_bjacprec type, extends(psb_s_base_prec_type) :: psb_s_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) type(psb_sspmat_type), allocatable :: av(:) - type(psb_s_vect_type), allocatable :: dv + type(psb_s_vect_type), allocatable :: dv, wrk(:) contains procedure, pass(prec) :: s_apply_v => psb_s_bjac_apply_vect procedure, pass(prec) :: s_apply => psb_s_bjac_apply @@ -49,6 +49,9 @@ module psb_s_bjacprec procedure, pass(prec) :: free => psb_s_bjac_precfree procedure, pass(prec) :: sizeof => psb_s_bjac_sizeof procedure, pass(prec) :: get_nzeros => psb_s_bjac_get_nzeros + procedure, pass(prec) :: allocate_wrk => psb_s_bjac_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_s_bjac_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_s_bjac_is_allocated_wrk end type psb_s_bjac_prec_type private :: psb_s_bjac_sizeof, psb_s_bjac_precdescr, psb_s_bjac_get_nzeros @@ -308,4 +311,111 @@ contains end subroutine psb_s_bjac_clone + subroutine psb_s_bjac_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_s_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act, i + character(len=20) :: name + + info=psb_success_ + name = 'psb_s_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + if (allocated(prec%wrk)) then + if (size(prec%wrk)<2) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + if (.not.allocated(prec%wrk)) then + if (.not.present(desc)) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 + end if + allocate(prec%wrk(2),stat=info) + do i=1, 2 + if (info == 0) call psb_geall(prec%wrk(i),desc,info) + if (info == 0) call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.) + end do + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_s_bjac_allocate_wrk + + subroutine psb_s_bjac_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_s_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: i + character(len=20) :: name + + info=psb_success_ + name = 'psb_s_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + info = psb_success_ + if (allocated(prec%wrk)) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_s_bjac_free_wrk + + function psb_s_bjac_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_s_bjac_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = allocated(prec%wrk) + + end function psb_s_bjac_is_allocated_wrk + + end module psb_s_bjacprec diff --git a/prec/psb_s_prec_type.f90 b/prec/psb_s_prec_type.f90 index 4eb157c6..ded50ca4 100644 --- a/prec/psb_s_prec_type.f90 +++ b/prec/psb_s_prec_type.f90 @@ -54,6 +54,9 @@ module psb_s_prec_type procedure, pass(prec) :: build => psb_sprecbld procedure, pass(prec) :: init => psb_sprecinit procedure, pass(prec) :: descr => psb_sfile_prec_descr + procedure, pass(prec) :: allocate_wrk => psb_s_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_s_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_s_is_allocated_wrk end type psb_sprec_type interface psb_precfree @@ -193,6 +196,90 @@ contains end subroutine psb_s_prec_dump + subroutine psb_s_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_sprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_s_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner' + return + end if + + call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_s_allocate_wrk + + subroutine psb_s_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_sprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_s_free_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to free a non-built preconditioner' + return + end if + + call prec%prec%free_wrk(info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_s_free_wrk + + function psb_s_is_allocated_wrk(prec) result(res) + implicit none + + ! Arguments + class(psb_sprec_type), intent(in) :: prec + logical :: res + + if (.not.allocated(prec%prec)) then + res = .false. + else + res = prec%prec%is_allocated_wrk() + end if + + end function psb_s_is_allocated_wrk subroutine psb_s_precfree(p,info) use psb_base_mod diff --git a/prec/psb_z_base_prec_mod.f90 b/prec/psb_z_base_prec_mod.f90 index ad9eacc1..a1e832f4 100644 --- a/prec/psb_z_base_prec_mod.f90 +++ b/prec/psb_z_base_prec_mod.f90 @@ -61,6 +61,9 @@ module psb_z_base_prec_mod generic, public :: build => precbld generic, public :: descr => precdescr procedure, pass(prec) :: desc_prefix => psb_z_base_desc_prefix + procedure, pass(prec) :: allocate_wrk => psb_z_base_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_z_base_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_z_base_is_allocated_wrk procedure(psb_z_base_precbld), pass(prec), deferred :: precbld procedure(psb_z_base_sizeof), pass(prec), deferred :: sizeof procedure(psb_z_base_precinit), pass(prec), deferred :: precinit @@ -258,6 +261,87 @@ contains end subroutine psb_z_base_precsetc + subroutine psb_z_base_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_z_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_z_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_z_base_allocate_wrk + + subroutine psb_z_base_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_z_base_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_z_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + ! + ! Base version does nothing. + ! + + info = psb_success_ + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_z_base_free_wrk + + function psb_z_base_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_z_base_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = .true. + + end function psb_z_base_is_allocated_wrk + subroutine psb_z_base_set_ctxt(prec,ictxt) implicit none class(psb_z_base_prec_type), intent(inout) :: prec diff --git a/prec/psb_z_bjacprec.f90 b/prec/psb_z_bjacprec.f90 index d7b630d0..8cad5fcd 100644 --- a/prec/psb_z_bjacprec.f90 +++ b/prec/psb_z_bjacprec.f90 @@ -36,7 +36,7 @@ module psb_z_bjacprec type, extends(psb_z_base_prec_type) :: psb_z_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) type(psb_zspmat_type), allocatable :: av(:) - type(psb_z_vect_type), allocatable :: dv + type(psb_z_vect_type), allocatable :: dv, wrk(:) contains procedure, pass(prec) :: z_apply_v => psb_z_bjac_apply_vect procedure, pass(prec) :: z_apply => psb_z_bjac_apply @@ -49,6 +49,9 @@ module psb_z_bjacprec procedure, pass(prec) :: free => psb_z_bjac_precfree procedure, pass(prec) :: sizeof => psb_z_bjac_sizeof procedure, pass(prec) :: get_nzeros => psb_z_bjac_get_nzeros + procedure, pass(prec) :: allocate_wrk => psb_z_bjac_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_z_bjac_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_z_bjac_is_allocated_wrk end type psb_z_bjac_prec_type private :: psb_z_bjac_sizeof, psb_z_bjac_precdescr, psb_z_bjac_get_nzeros @@ -308,4 +311,111 @@ contains end subroutine psb_z_bjac_clone + subroutine psb_z_bjac_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_z_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act, i + character(len=20) :: name + + info=psb_success_ + name = 'psb_z_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + if (allocated(prec%wrk)) then + if (size(prec%wrk)<2) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + if (.not.allocated(prec%wrk)) then + if (.not.present(desc)) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 + end if + allocate(prec%wrk(2),stat=info) + do i=1, 2 + if (info == 0) call psb_geall(prec%wrk(i),desc,info) + if (info == 0) call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.) + end do + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_z_bjac_allocate_wrk + + subroutine psb_z_bjac_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_z_bjac_prec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: i + character(len=20) :: name + + info=psb_success_ + name = 'psb_z_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + info = psb_success_ + if (allocated(prec%wrk)) then + do i=1,size(prec%wrk) + if (info == 0) call prec%wrk(i)%free(info) + end do + if (info == 0) deallocate(prec%wrk,stat=info) + end if + if (info /= 0) then + info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_z_bjac_free_wrk + + function psb_z_bjac_is_allocated_wrk(prec) result(res) + use psb_base_mod + implicit none + + ! Arguments + class(psb_z_bjac_prec_type), intent(in) :: prec + logical :: res + + ! In the base version we can say yes, because + ! there is nothing to allocate + + res = allocated(prec%wrk) + + end function psb_z_bjac_is_allocated_wrk + + end module psb_z_bjacprec diff --git a/prec/psb_z_prec_type.f90 b/prec/psb_z_prec_type.f90 index 15a53fe4..8b1a8f4b 100644 --- a/prec/psb_z_prec_type.f90 +++ b/prec/psb_z_prec_type.f90 @@ -54,6 +54,9 @@ module psb_z_prec_type procedure, pass(prec) :: build => psb_zprecbld procedure, pass(prec) :: init => psb_zprecinit procedure, pass(prec) :: descr => psb_zfile_prec_descr + procedure, pass(prec) :: allocate_wrk => psb_z_allocate_wrk + procedure, pass(prec) :: free_wrk => psb_z_free_wrk + procedure, pass(prec) :: is_allocated_wrk => psb_z_is_allocated_wrk end type psb_zprec_type interface psb_precfree @@ -193,6 +196,90 @@ contains end subroutine psb_z_prec_dump + subroutine psb_z_allocate_wrk(prec,info,vmold,desc) + use psb_base_mod + implicit none + + ! Arguments + class(psb_zprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_vect_type), intent(in), optional :: vmold + type(psb_desc_type), intent(in), optional :: desc + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_z_allocate_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner' + return + end if + + call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_z_allocate_wrk + + subroutine psb_z_free_wrk(prec,info) + use psb_base_mod + implicit none + + ! Arguments + class(psb_zprec_type), intent(inout) :: prec + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20) :: name + + info=psb_success_ + name = 'psb_z_free_wrk' + call psb_erractionsave(err_act) + + if (psb_get_errstatus().ne.0) goto 9999 + + if (.not.allocated(prec%prec)) then + info = -1 + write(psb_err_unit,*) 'Trying to free a non-built preconditioner' + return + end if + + call prec%prec%free_wrk(info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine psb_z_free_wrk + + function psb_z_is_allocated_wrk(prec) result(res) + implicit none + + ! Arguments + class(psb_zprec_type), intent(in) :: prec + logical :: res + + if (.not.allocated(prec%prec)) then + res = .false. + else + res = prec%prec%is_allocated_wrk() + end if + + end function psb_z_is_allocated_wrk subroutine psb_z_precfree(p,info) use psb_base_mod