From 62c75abbf493c2bea21732bf377268d2036db609 Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Mon, 16 Nov 2020 18:37:45 +0100 Subject: [PATCH] Completed integration of ILU-type factorization --- prec/impl/psb_c_bjacprec_impl.f90 | 44 +---- prec/impl/psb_c_prec_type_impl.f90 | 225 +++++++++++++++++++------ prec/impl/psb_d_bjacprec_impl.f90 | 44 +---- prec/impl/psb_d_prec_type_impl.f90 | 225 +++++++++++++++++++------ prec/impl/psb_s_bjacprec_impl.f90 | 44 +---- prec/impl/psb_s_prec_type_impl.f90 | 225 +++++++++++++++++++------ prec/impl/psb_z_bjacprec_impl.f90 | 44 +---- prec/impl/psb_z_prec_type_impl.f90 | 225 +++++++++++++++++++------ prec/psb_c_bjacprec.f90 | 8 +- prec/psb_c_prec_type.f90 | 151 ++++++++++------- prec/psb_d_bjacprec.f90 | 8 +- prec/psb_d_prec_type.f90 | 151 ++++++++++------- prec/psb_s_bjacprec.f90 | 8 +- prec/psb_s_prec_type.f90 | 151 ++++++++++------- prec/psb_z_bjacprec.f90 | 8 +- prec/psb_z_prec_type.f90 | 151 ++++++++++------- test/pargen/psb_d_pde2d.f90 | 254 ++++++++++++++++++----------- test/pargen/psb_d_pde3d.f90 | 61 +++++-- test/pargen/psb_s_pde2d.f90 | 254 ++++++++++++++++++----------- test/pargen/psb_s_pde3d.f90 | 63 +++++-- test/pargen/runs/ppde.inp | 18 +- 21 files changed, 1510 insertions(+), 852 deletions(-) diff --git a/prec/impl/psb_c_bjacprec_impl.f90 b/prec/impl/psb_c_bjacprec_impl.f90 index 3cef79dc..02e23da3 100644 --- a/prec/impl/psb_c_bjacprec_impl.f90 +++ b/prec/impl/psb_c_bjacprec_impl.f90 @@ -568,7 +568,6 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 endif ! This is where we have no renumbering, thus no need - ! call psb_ilu_fct(a,lf,uf,dd,info) call psb_ilu0_fact(ialg,a,lf,uf,dd,info) if(info == psb_success_) then @@ -782,45 +781,19 @@ subroutine psb_c_bjac_precseti(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_ilu_fill_in_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_fill_in_) = val case (psb_ilu_ialg_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_ialg_) = val case (psb_ilu_scale_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_scale_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select @@ -855,26 +828,13 @@ subroutine psb_c_bjac_precsetr(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_fact_eps_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%rprcparm(psb_fact_eps_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select diff --git a/prec/impl/psb_c_prec_type_impl.f90 b/prec/impl/psb_c_prec_type_impl.f90 index e1f13fc4..1d99a85d 100644 --- a/prec/impl/psb_c_prec_type_impl.f90 +++ b/prec/impl/psb_c_prec_type_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,14 +27,14 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! -! +! +! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -46,7 +46,7 @@ !!$ 3. The name of the PSBLAS group or the names of its contributors may !!$ not be used to endorse or promote products derived from this !!$ software without specific written permission. -!!$ +!!$ !!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS !!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED !!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -58,14 +58,14 @@ !!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +!!$ +!!$ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine psb_c_apply2_vect(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_c_prec_type, psb_protect_name => psb_c_apply2_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_cprec_type), intent(inout) :: prec type(psb_c_vect_type),intent(inout) :: x @@ -74,7 +74,7 @@ subroutine psb_c_apply2_vect(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans complex(psb_spk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ complex(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -87,25 +87,25 @@ subroutine psb_c_apply2_vect(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -114,13 +114,13 @@ subroutine psb_c_apply2_vect(prec,x,y,desc_data,info,trans,work) call prec%prec%apply(cone,x,czero,y,desc_data,info,& & trans=trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -135,7 +135,7 @@ end subroutine psb_c_apply2_vect subroutine psb_c_apply1_vect(prec,x,desc_data,info,trans,work) use psb_base_mod use psb_c_prec_type, psb_protect_name => psb_c_apply1_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_cprec_type), intent(inout) :: prec type(psb_c_vect_type),intent(inout) :: x @@ -144,7 +144,7 @@ subroutine psb_c_apply1_vect(prec,x,desc_data,info,trans,work) complex(psb_spk_),intent(inout), optional, target :: work(:) type(psb_c_vect_type) :: ww - character :: trans_ + character :: trans_ complex(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -157,25 +157,25 @@ subroutine psb_c_apply1_vect(prec,x,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -186,13 +186,13 @@ subroutine psb_c_apply1_vect(prec,x,desc_data,info,trans,work) & trans=trans_,work=work_) if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info) call psb_gefree(ww,desc_data,info) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -207,7 +207,7 @@ end subroutine psb_c_apply1_vect subroutine psb_c_apply2v(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_c_prec_type, psb_protect_name => psb_c_apply2v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_cprec_type), intent(inout) :: prec complex(psb_spk_),intent(inout) :: x(:) @@ -216,7 +216,7 @@ subroutine psb_c_apply2v(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans complex(psb_spk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ complex(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -229,37 +229,37 @@ subroutine psb_c_apply2v(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=trans else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if call prec%prec%apply(cone,x,czero,y,desc_data,info,trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -274,7 +274,7 @@ end subroutine psb_c_apply2v subroutine psb_c_apply1v(prec,x,desc_data,info,trans) use psb_base_mod use psb_c_prec_type, psb_protect_name => psb_c_apply1v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_cprec_type), intent(inout) :: prec complex(psb_spk_),intent(inout) :: x(:) @@ -293,32 +293,32 @@ subroutine psb_c_apply1v(prec,x,desc_data,info,trans) ictxt=desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if call prec%prec%apply(cone,x,czero,ww,desc_data,info,& & trans_,work=w1) if(info /= psb_success_) goto 9999 x(:) = ww(:) deallocate(ww,W1,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if @@ -332,3 +332,126 @@ subroutine psb_c_apply1v(prec,x,desc_data,info,trans) end subroutine psb_c_apply1v +subroutine psb_ccprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_c_prec_type, psb_protect_name => psb_ccprecseti + implicit none + + class(psb_cprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='psb_precseti' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ("SUB_FILLIN") + call prec%prec%precset(psb_ilu_fill_in_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_ccprecseti + +subroutine psb_ccprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_c_prec_type, psb_protect_name => psb_ccprecsetr + implicit none + + class(psb_cprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetr' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case('SUB_ILUTHRS') + call prec%prec%precset(psb_fact_eps_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_ccprecsetr + +subroutine psb_ccprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_c_prec_type, psb_protect_name => psb_ccprecsetc + implicit none + + class(psb_cprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetc' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ('SUB_SOLVE') + ! We select here the type of solver on the block + select case (psb_toupper(string)) + case("ILU") + call prec%prec%precset(psb_f_type_,psb_f_ilu_k_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + case("ILUT") + call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case default + ! Default to ILU(0) factorization + call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + end select + case ("ILU_ALG") + select case (psb_toupper(string)) + case ("MILU") + call prec%prec%precset(psb_ilu_ialg_,psb_milu_n_,info) + case default + ! Do nothing + end select + case ("ILUT_SCALE") + select case (psb_toupper(string)) + case ("MAXVAL") + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info) + case default + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) + end select + case default + + end select + +end subroutine psb_ccprecsetc diff --git a/prec/impl/psb_d_bjacprec_impl.f90 b/prec/impl/psb_d_bjacprec_impl.f90 index 34a2b63a..898fd224 100644 --- a/prec/impl/psb_d_bjacprec_impl.f90 +++ b/prec/impl/psb_d_bjacprec_impl.f90 @@ -568,7 +568,6 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 endif ! This is where we have no renumbering, thus no need - ! call psb_ilu_fct(a,lf,uf,dd,info) call psb_ilu0_fact(ialg,a,lf,uf,dd,info) if(info == psb_success_) then @@ -782,45 +781,19 @@ subroutine psb_d_bjac_precseti(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_ilu_fill_in_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_fill_in_) = val case (psb_ilu_ialg_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_ialg_) = val case (psb_ilu_scale_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_scale_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select @@ -855,26 +828,13 @@ subroutine psb_d_bjac_precsetr(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_fact_eps_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%rprcparm(psb_fact_eps_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select diff --git a/prec/impl/psb_d_prec_type_impl.f90 b/prec/impl/psb_d_prec_type_impl.f90 index 793afac7..5c07bebd 100644 --- a/prec/impl/psb_d_prec_type_impl.f90 +++ b/prec/impl/psb_d_prec_type_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,14 +27,14 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! -! +! +! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -46,7 +46,7 @@ !!$ 3. The name of the PSBLAS group or the names of its contributors may !!$ not be used to endorse or promote products derived from this !!$ software without specific written permission. -!!$ +!!$ !!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS !!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED !!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -58,14 +58,14 @@ !!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +!!$ +!!$ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine psb_d_apply2_vect(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_d_prec_type, psb_protect_name => psb_d_apply2_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_dprec_type), intent(inout) :: prec type(psb_d_vect_type),intent(inout) :: x @@ -74,7 +74,7 @@ subroutine psb_d_apply2_vect(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans real(psb_dpk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ real(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -87,25 +87,25 @@ subroutine psb_d_apply2_vect(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -114,13 +114,13 @@ subroutine psb_d_apply2_vect(prec,x,y,desc_data,info,trans,work) call prec%prec%apply(done,x,dzero,y,desc_data,info,& & trans=trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -135,7 +135,7 @@ end subroutine psb_d_apply2_vect subroutine psb_d_apply1_vect(prec,x,desc_data,info,trans,work) use psb_base_mod use psb_d_prec_type, psb_protect_name => psb_d_apply1_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_dprec_type), intent(inout) :: prec type(psb_d_vect_type),intent(inout) :: x @@ -144,7 +144,7 @@ subroutine psb_d_apply1_vect(prec,x,desc_data,info,trans,work) real(psb_dpk_),intent(inout), optional, target :: work(:) type(psb_d_vect_type) :: ww - character :: trans_ + character :: trans_ real(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -157,25 +157,25 @@ subroutine psb_d_apply1_vect(prec,x,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -186,13 +186,13 @@ subroutine psb_d_apply1_vect(prec,x,desc_data,info,trans,work) & trans=trans_,work=work_) if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info) call psb_gefree(ww,desc_data,info) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -207,7 +207,7 @@ end subroutine psb_d_apply1_vect subroutine psb_d_apply2v(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_d_prec_type, psb_protect_name => psb_d_apply2v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_dprec_type), intent(inout) :: prec real(psb_dpk_),intent(inout) :: x(:) @@ -216,7 +216,7 @@ subroutine psb_d_apply2v(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans real(psb_dpk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ real(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -229,37 +229,37 @@ subroutine psb_d_apply2v(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=trans else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if call prec%prec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -274,7 +274,7 @@ end subroutine psb_d_apply2v subroutine psb_d_apply1v(prec,x,desc_data,info,trans) use psb_base_mod use psb_d_prec_type, psb_protect_name => psb_d_apply1v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_dprec_type), intent(inout) :: prec real(psb_dpk_),intent(inout) :: x(:) @@ -293,32 +293,32 @@ subroutine psb_d_apply1v(prec,x,desc_data,info,trans) ictxt=desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if call prec%prec%apply(done,x,dzero,ww,desc_data,info,& & trans_,work=w1) if(info /= psb_success_) goto 9999 x(:) = ww(:) deallocate(ww,W1,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if @@ -332,3 +332,126 @@ subroutine psb_d_apply1v(prec,x,desc_data,info,trans) end subroutine psb_d_apply1v +subroutine psb_dcprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_d_prec_type, psb_protect_name => psb_dcprecseti + implicit none + + class(psb_dprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='psb_precseti' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ("SUB_FILLIN") + call prec%prec%precset(psb_ilu_fill_in_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_dcprecseti + +subroutine psb_dcprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_d_prec_type, psb_protect_name => psb_dcprecsetr + implicit none + + class(psb_dprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetr' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case('SUB_ILUTHRS') + call prec%prec%precset(psb_fact_eps_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_dcprecsetr + +subroutine psb_dcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_d_prec_type, psb_protect_name => psb_dcprecsetc + implicit none + + class(psb_dprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetc' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ('SUB_SOLVE') + ! We select here the type of solver on the block + select case (psb_toupper(string)) + case("ILU") + call prec%prec%precset(psb_f_type_,psb_f_ilu_k_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + case("ILUT") + call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case default + ! Default to ILU(0) factorization + call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + end select + case ("ILU_ALG") + select case (psb_toupper(string)) + case ("MILU") + call prec%prec%precset(psb_ilu_ialg_,psb_milu_n_,info) + case default + ! Do nothing + end select + case ("ILUT_SCALE") + select case (psb_toupper(string)) + case ("MAXVAL") + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info) + case default + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) + end select + case default + + end select + +end subroutine psb_dcprecsetc diff --git a/prec/impl/psb_s_bjacprec_impl.f90 b/prec/impl/psb_s_bjacprec_impl.f90 index 1448e43f..b2545890 100644 --- a/prec/impl/psb_s_bjacprec_impl.f90 +++ b/prec/impl/psb_s_bjacprec_impl.f90 @@ -568,7 +568,6 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 endif ! This is where we have no renumbering, thus no need - ! call psb_ilu_fct(a,lf,uf,dd,info) call psb_ilu0_fact(ialg,a,lf,uf,dd,info) if(info == psb_success_) then @@ -782,45 +781,19 @@ subroutine psb_s_bjac_precseti(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_ilu_fill_in_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_fill_in_) = val case (psb_ilu_ialg_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_ialg_) = val case (psb_ilu_scale_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_scale_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select @@ -855,26 +828,13 @@ subroutine psb_s_bjac_precsetr(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_fact_eps_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%rprcparm(psb_fact_eps_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select diff --git a/prec/impl/psb_s_prec_type_impl.f90 b/prec/impl/psb_s_prec_type_impl.f90 index 547272a0..837c3106 100644 --- a/prec/impl/psb_s_prec_type_impl.f90 +++ b/prec/impl/psb_s_prec_type_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,14 +27,14 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! -! +! +! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -46,7 +46,7 @@ !!$ 3. The name of the PSBLAS group or the names of its contributors may !!$ not be used to endorse or promote products derived from this !!$ software without specific written permission. -!!$ +!!$ !!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS !!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED !!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -58,14 +58,14 @@ !!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +!!$ +!!$ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine psb_s_apply2_vect(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_s_prec_type, psb_protect_name => psb_s_apply2_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_sprec_type), intent(inout) :: prec type(psb_s_vect_type),intent(inout) :: x @@ -74,7 +74,7 @@ subroutine psb_s_apply2_vect(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans real(psb_spk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ real(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -87,25 +87,25 @@ subroutine psb_s_apply2_vect(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -114,13 +114,13 @@ subroutine psb_s_apply2_vect(prec,x,y,desc_data,info,trans,work) call prec%prec%apply(sone,x,szero,y,desc_data,info,& & trans=trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -135,7 +135,7 @@ end subroutine psb_s_apply2_vect subroutine psb_s_apply1_vect(prec,x,desc_data,info,trans,work) use psb_base_mod use psb_s_prec_type, psb_protect_name => psb_s_apply1_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_sprec_type), intent(inout) :: prec type(psb_s_vect_type),intent(inout) :: x @@ -144,7 +144,7 @@ subroutine psb_s_apply1_vect(prec,x,desc_data,info,trans,work) real(psb_spk_),intent(inout), optional, target :: work(:) type(psb_s_vect_type) :: ww - character :: trans_ + character :: trans_ real(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -157,25 +157,25 @@ subroutine psb_s_apply1_vect(prec,x,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -186,13 +186,13 @@ subroutine psb_s_apply1_vect(prec,x,desc_data,info,trans,work) & trans=trans_,work=work_) if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info) call psb_gefree(ww,desc_data,info) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -207,7 +207,7 @@ end subroutine psb_s_apply1_vect subroutine psb_s_apply2v(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_s_prec_type, psb_protect_name => psb_s_apply2v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_sprec_type), intent(inout) :: prec real(psb_spk_),intent(inout) :: x(:) @@ -216,7 +216,7 @@ subroutine psb_s_apply2v(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans real(psb_spk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ real(psb_spk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -229,37 +229,37 @@ subroutine psb_s_apply2v(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=trans else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if call prec%prec%apply(sone,x,szero,y,desc_data,info,trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -274,7 +274,7 @@ end subroutine psb_s_apply2v subroutine psb_s_apply1v(prec,x,desc_data,info,trans) use psb_base_mod use psb_s_prec_type, psb_protect_name => psb_s_apply1v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_sprec_type), intent(inout) :: prec real(psb_spk_),intent(inout) :: x(:) @@ -293,32 +293,32 @@ subroutine psb_s_apply1v(prec,x,desc_data,info,trans) ictxt=desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if call prec%prec%apply(sone,x,szero,ww,desc_data,info,& & trans_,work=w1) if(info /= psb_success_) goto 9999 x(:) = ww(:) deallocate(ww,W1,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if @@ -332,3 +332,126 @@ subroutine psb_s_apply1v(prec,x,desc_data,info,trans) end subroutine psb_s_apply1v +subroutine psb_scprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_s_prec_type, psb_protect_name => psb_scprecseti + implicit none + + class(psb_sprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='psb_precseti' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ("SUB_FILLIN") + call prec%prec%precset(psb_ilu_fill_in_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_scprecseti + +subroutine psb_scprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_s_prec_type, psb_protect_name => psb_scprecsetr + implicit none + + class(psb_sprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetr' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case('SUB_ILUTHRS') + call prec%prec%precset(psb_fact_eps_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_scprecsetr + +subroutine psb_scprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_s_prec_type, psb_protect_name => psb_scprecsetc + implicit none + + class(psb_sprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetc' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ('SUB_SOLVE') + ! We select here the type of solver on the block + select case (psb_toupper(string)) + case("ILU") + call prec%prec%precset(psb_f_type_,psb_f_ilu_k_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + case("ILUT") + call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case default + ! Default to ILU(0) factorization + call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + end select + case ("ILU_ALG") + select case (psb_toupper(string)) + case ("MILU") + call prec%prec%precset(psb_ilu_ialg_,psb_milu_n_,info) + case default + ! Do nothing + end select + case ("ILUT_SCALE") + select case (psb_toupper(string)) + case ("MAXVAL") + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info) + case default + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) + end select + case default + + end select + +end subroutine psb_scprecsetc diff --git a/prec/impl/psb_z_bjacprec_impl.f90 b/prec/impl/psb_z_bjacprec_impl.f90 index a1859e87..2a8960ea 100644 --- a/prec/impl/psb_z_bjacprec_impl.f90 +++ b/prec/impl/psb_z_bjacprec_impl.f90 @@ -568,7 +568,6 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 endif ! This is where we have no renumbering, thus no need - ! call psb_ilu_fct(a,lf,uf,dd,info) call psb_ilu0_fact(ialg,a,lf,uf,dd,info) if(info == psb_success_) then @@ -782,45 +781,19 @@ subroutine psb_z_bjac_precseti(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_ilu_fill_in_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_fill_in_) = val case (psb_ilu_ialg_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_ialg_) = val case (psb_ilu_scale_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_ilu_scale_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select @@ -855,26 +828,13 @@ subroutine psb_z_bjac_precsetr(prec,what,val,info) select case(what) case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%iprcparm(psb_f_type_) = val case (psb_fact_eps_) - if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then - write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& - & prec%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif prec%rprcparm(psb_fact_eps_) = val case default - write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + write(psb_err_unit,'(i0," is invalid, ignoring user specification")') what end select diff --git a/prec/impl/psb_z_prec_type_impl.f90 b/prec/impl/psb_z_prec_type_impl.f90 index 982fc008..3f4afe66 100644 --- a/prec/impl/psb_z_prec_type_impl.f90 +++ b/prec/impl/psb_z_prec_type_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,14 +27,14 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! -! +! +! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -46,7 +46,7 @@ !!$ 3. The name of the PSBLAS group or the names of its contributors may !!$ not be used to endorse or promote products derived from this !!$ software without specific written permission. -!!$ +!!$ !!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS !!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED !!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -58,14 +58,14 @@ !!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) !!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE !!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ +!!$ +!!$ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine psb_z_apply2_vect(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_z_prec_type, psb_protect_name => psb_z_apply2_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_zprec_type), intent(inout) :: prec type(psb_z_vect_type),intent(inout) :: x @@ -74,7 +74,7 @@ subroutine psb_z_apply2_vect(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans complex(psb_dpk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ complex(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -87,25 +87,25 @@ subroutine psb_z_apply2_vect(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -114,13 +114,13 @@ subroutine psb_z_apply2_vect(prec,x,y,desc_data,info,trans,work) call prec%prec%apply(zone,x,zzero,y,desc_data,info,& & trans=trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -135,7 +135,7 @@ end subroutine psb_z_apply2_vect subroutine psb_z_apply1_vect(prec,x,desc_data,info,trans,work) use psb_base_mod use psb_z_prec_type, psb_protect_name => psb_z_apply1_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_zprec_type), intent(inout) :: prec type(psb_z_vect_type),intent(inout) :: x @@ -144,7 +144,7 @@ subroutine psb_z_apply1_vect(prec,x,desc_data,info,trans,work) complex(psb_dpk_),intent(inout), optional, target :: work(:) type(psb_z_vect_type) :: ww - character :: trans_ + character :: trans_ complex(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -157,25 +157,25 @@ subroutine psb_z_apply1_vect(prec,x,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -186,13 +186,13 @@ subroutine psb_z_apply1_vect(prec,x,desc_data,info,trans,work) & trans=trans_,work=work_) if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info) call psb_gefree(ww,desc_data,info) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -207,7 +207,7 @@ end subroutine psb_z_apply1_vect subroutine psb_z_apply2v(prec,x,y,desc_data,info,trans,work) use psb_base_mod use psb_z_prec_type, psb_protect_name => psb_z_apply2v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_zprec_type), intent(inout) :: prec complex(psb_dpk_),intent(inout) :: x(:) @@ -216,7 +216,7 @@ subroutine psb_z_apply2v(prec,x,y,desc_data,info,trans,work) character(len=1), optional :: trans complex(psb_dpk_),intent(inout), optional, target :: work(:) - character :: trans_ + character :: trans_ complex(psb_dpk_), pointer :: work_(:) integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: err_act @@ -229,37 +229,37 @@ subroutine psb_z_apply2v(prec,x,y,desc_data,info,trans,work) ictxt = desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=trans else trans_='N' end if - if (present(work)) then + if (present(work)) then work_ => work else allocate(work_(4*desc_data%get_local_cols()),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if call prec%prec%apply(zone,x,zzero,y,desc_data,info,trans_,work=work_) - if (present(work)) then + if (present(work)) then else deallocate(work_,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if end if @@ -274,7 +274,7 @@ end subroutine psb_z_apply2v subroutine psb_z_apply1v(prec,x,desc_data,info,trans) use psb_base_mod use psb_z_prec_type, psb_protect_name => psb_z_apply1v - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_zprec_type), intent(inout) :: prec complex(psb_dpk_),intent(inout) :: x(:) @@ -293,32 +293,32 @@ subroutine psb_z_apply1v(prec,x,desc_data,info,trans) ictxt=desc_data%get_context() call psb_info(ictxt, me, np) - if (present(trans)) then + if (present(trans)) then trans_=psb_toupper(trans) else trans_='N' end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') - goto 9999 + goto 9999 end if call prec%prec%apply(zone,x,zzero,ww,desc_data,info,& & trans_,work=w1) if(info /= psb_success_) goto 9999 x(:) = ww(:) deallocate(ww,W1,stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then info = psb_err_from_subroutine_ call psb_errpush(info,name,a_err='DeAllocate') - goto 9999 + goto 9999 end if @@ -332,3 +332,126 @@ subroutine psb_z_apply1v(prec,x,desc_data,info,trans) end subroutine psb_z_apply1v +subroutine psb_zcprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_z_prec_type, psb_protect_name => psb_zcprecseti + implicit none + + class(psb_zprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='psb_precseti' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ("SUB_FILLIN") + call prec%prec%precset(psb_ilu_fill_in_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_zcprecseti + +subroutine psb_zcprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_z_prec_type, psb_protect_name => psb_zcprecsetr + implicit none + + class(psb_zprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetr' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case('SUB_ILUTHRS') + call prec%prec%precset(psb_fact_eps_,val,info) + case default + info = psb_err_invalid_args_combination_ + write(psb_err_unit,*) name,& + & ': Error: uninitialized preconditioner,',& + &' should call prec%init' + return + end select + +end subroutine psb_zcprecsetr + +subroutine psb_zcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + use psb_base_mod + use psb_z_prec_type, psb_protect_name => psb_zcprecsetc + implicit none + + class(psb_zprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + ! This optional inputs are backport from the inputs available in AMG4PSBLAS, + ! they are of no actual use here a part from compatibility reasons. + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + + ! Local variables + character(len=*), parameter :: name='amg_precsetc' + + info = psb_success_ + + ! We need to convert from the 'what' string to the corresponding integer + ! value befor passing the call to the set of the inner method. + select case (psb_toupper(what)) + case ('SUB_SOLVE') + ! We select here the type of solver on the block + select case (psb_toupper(string)) + case("ILU") + call prec%prec%precset(psb_f_type_,psb_f_ilu_k_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + case("ILUT") + call prec%prec%precset(psb_f_type_,psb_f_ilu_t_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_t_,info) + case default + ! Default to ILU(0) factorization + call prec%prec%precset(psb_f_type_,psb_f_ilu_n_,info) + call prec%prec%precset(psb_ilu_ialg_,psb_ilu_n_,info) + end select + case ("ILU_ALG") + select case (psb_toupper(string)) + case ("MILU") + call prec%prec%precset(psb_ilu_ialg_,psb_milu_n_,info) + case default + ! Do nothing + end select + case ("ILUT_SCALE") + select case (psb_toupper(string)) + case ("MAXVAL") + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info) + case default + call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info) + end select + case default + + end select + +end subroutine psb_zcprecsetc diff --git a/prec/psb_c_bjacprec.f90 b/prec/psb_c_bjacprec.f90 index 7a232c46..ca4cbcb8 100644 --- a/prec/psb_c_bjacprec.f90 +++ b/prec/psb_c_bjacprec.f90 @@ -57,13 +57,11 @@ module psb_c_bjacprec 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 - - character(len=15), parameter, private :: & - & fact_names(0:2)=(/'None ','ILU(n) ',& - & 'ILU(eps) '/) + & fact_names(0:3)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) '/) + private :: psb_c_bjac_sizeof, psb_c_bjac_precdescr, psb_c_bjac_get_nzeros interface subroutine psb_c_bjac_dump(prec,info,prefix,head) diff --git a/prec/psb_c_prec_type.f90 b/prec/psb_c_prec_type.f90 index 6d9aa908..ce03f635 100644 --- a/prec/psb_c_prec_type.f90 +++ b/prec/psb_c_prec_type.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,8 +27,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Module to define PREC_DATA, !! !! structure for preconditioning. !! @@ -39,7 +39,7 @@ module psb_c_prec_type use psb_c_base_prec_mod type psb_cprec_type - integer(psb_ipk_) :: ictxt + integer(psb_ipk_) :: ictxt class(psb_c_base_prec_type), allocatable :: prec contains procedure, pass(prec) :: psb_c_apply1_vect @@ -54,6 +54,10 @@ 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) :: cseti => psb_ccprecseti + procedure, pass(prec) :: csetc => psb_ccprecsetc + procedure, pass(prec) :: csetr => psb_ccprecsetr + generic, public :: set => cseti, csetc, csetr 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 @@ -102,7 +106,7 @@ module psb_c_prec_type module procedure psb_cprec_sizeof end interface - interface + interface subroutine psb_c_apply2_vect(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_cprec_type, psb_c_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -114,8 +118,8 @@ module psb_c_prec_type complex(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_c_apply2_vect end interface - - interface + + interface subroutine psb_c_apply1_vect(prec,x,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_cprec_type, psb_c_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -126,7 +130,7 @@ module psb_c_prec_type complex(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_c_apply1_vect end interface - + interface subroutine psb_c_apply2v(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_cprec_type, psb_c_vect_type, psb_spk_ @@ -139,8 +143,8 @@ module psb_c_prec_type complex(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_c_apply2v end interface - - interface + + interface subroutine psb_c_apply1v(prec,x,desc_data,info,trans) import :: psb_ipk_, psb_desc_type, psb_cprec_type, psb_c_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -150,56 +154,89 @@ module psb_c_prec_type character(len=1), optional :: trans end subroutine psb_c_apply1v end interface - + + interface + subroutine psb_ccprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_cprec_type, psb_cspmat_type, psb_desc_type, psb_spk_, & + & psb_ipk_ + class(psb_cprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_ccprecseti + subroutine psb_ccprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_cprec_type, psb_cspmat_type, psb_desc_type, psb_spk_, & + & psb_ipk_ + class(psb_cprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_ccprecsetr + subroutine psb_ccprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + import :: psb_cprec_type, psb_cspmat_type, psb_desc_type, psb_spk_, & + & psb_ipk_ + class(psb_cprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_ccprecsetc +end interface + contains subroutine psb_cfile_prec_descr(prec,iout, root) use psb_base_mod - implicit none + implicit none class(psb_cprec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: root integer(psb_ipk_) :: iout_,info - character(len=20) :: name='prec_descr' - - if (present(iout)) then + character(len=20) :: name='prec_descr' + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") end if call prec%prec%descr(iout=iout,root=root) - + end subroutine psb_cfile_prec_descr subroutine psb_c_prec_dump(prec,info,prefix,head) - implicit none + implicit none type(psb_cprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head - ! len of prefix_ + ! len of prefix_ info = 0 - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = -1 write(psb_err_unit,*) 'Trying to dump a non-built preconditioner' return end if - + call prec%prec%dump(info,prefix,head) - - + + 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 @@ -209,33 +246,33 @@ contains ! 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 + 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 @@ -243,47 +280,47 @@ contains ! 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 + 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 + 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 - implicit none + implicit none type(psb_cprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -303,14 +340,14 @@ contains return 9999 call psb_error_handler(err_act) - + return end subroutine psb_c_precfree subroutine psb_c_prec_free(prec,info) use psb_base_mod - implicit none + implicit none class(psb_cprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -324,7 +361,7 @@ contains me=-1 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%free(info) if (info /= psb_success_) goto 9999 deallocate(prec%prec,stat=info) @@ -343,26 +380,26 @@ contains class(psb_cprec_type), intent(in) :: prec integer(psb_epk_) :: val integer(psb_ipk_) :: i - + val = 0 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then val = val + prec%prec%sizeof() end if - + end function psb_cprec_sizeof subroutine psb_c_prec_clone(prec,precout,info) - implicit none + implicit none class(psb_cprec_type), intent(inout) :: prec class(psb_cprec_type), intent(inout) :: precout integer(psb_ipk_), intent(out) :: info info = psb_success_ call prec%free(info) - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%clone(precout%prec,info) end if - + end subroutine psb_c_prec_clone end module psb_c_prec_type diff --git a/prec/psb_d_bjacprec.f90 b/prec/psb_d_bjacprec.f90 index e8c8e47a..0452d148 100644 --- a/prec/psb_d_bjacprec.f90 +++ b/prec/psb_d_bjacprec.f90 @@ -57,13 +57,11 @@ module psb_d_bjacprec 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 - - character(len=15), parameter, private :: & - & fact_names(0:2)=(/'None ','ILU(n) ',& - & 'ILU(eps) '/) + & fact_names(0:3)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) '/) + private :: psb_d_bjac_sizeof, psb_d_bjac_precdescr, psb_d_bjac_get_nzeros interface subroutine psb_d_bjac_dump(prec,info,prefix,head) diff --git a/prec/psb_d_prec_type.f90 b/prec/psb_d_prec_type.f90 index f50339a3..1b5b0f34 100644 --- a/prec/psb_d_prec_type.f90 +++ b/prec/psb_d_prec_type.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,8 +27,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Module to define PREC_DATA, !! !! structure for preconditioning. !! @@ -39,7 +39,7 @@ module psb_d_prec_type use psb_d_base_prec_mod type psb_dprec_type - integer(psb_ipk_) :: ictxt + integer(psb_ipk_) :: ictxt class(psb_d_base_prec_type), allocatable :: prec contains procedure, pass(prec) :: psb_d_apply1_vect @@ -54,6 +54,10 @@ 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) :: cseti => psb_dcprecseti + procedure, pass(prec) :: csetc => psb_dcprecsetc + procedure, pass(prec) :: csetr => psb_dcprecsetr + generic, public :: set => cseti, csetc, csetr 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 @@ -102,7 +106,7 @@ module psb_d_prec_type module procedure psb_dprec_sizeof end interface - interface + interface subroutine psb_d_apply2_vect(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_dprec_type, psb_d_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -114,8 +118,8 @@ module psb_d_prec_type real(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_d_apply2_vect end interface - - interface + + interface subroutine psb_d_apply1_vect(prec,x,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_dprec_type, psb_d_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -126,7 +130,7 @@ module psb_d_prec_type real(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_d_apply1_vect end interface - + interface subroutine psb_d_apply2v(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_dprec_type, psb_d_vect_type, psb_dpk_ @@ -139,8 +143,8 @@ module psb_d_prec_type real(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_d_apply2v end interface - - interface + + interface subroutine psb_d_apply1v(prec,x,desc_data,info,trans) import :: psb_ipk_, psb_desc_type, psb_dprec_type, psb_d_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -150,56 +154,89 @@ module psb_d_prec_type character(len=1), optional :: trans end subroutine psb_d_apply1v end interface - + + interface + subroutine psb_dcprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_dprec_type, psb_dspmat_type, psb_desc_type, psb_dpk_, & + & psb_ipk_ + class(psb_dprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_dcprecseti + subroutine psb_dcprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_dprec_type, psb_dspmat_type, psb_desc_type, psb_dpk_, & + & psb_ipk_ + class(psb_dprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_dcprecsetr + subroutine psb_dcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + import :: psb_dprec_type, psb_dspmat_type, psb_desc_type, psb_dpk_, & + & psb_ipk_ + class(psb_dprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_dcprecsetc +end interface + contains subroutine psb_dfile_prec_descr(prec,iout, root) use psb_base_mod - implicit none + implicit none class(psb_dprec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: root integer(psb_ipk_) :: iout_,info - character(len=20) :: name='prec_descr' - - if (present(iout)) then + character(len=20) :: name='prec_descr' + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") end if call prec%prec%descr(iout=iout,root=root) - + end subroutine psb_dfile_prec_descr subroutine psb_d_prec_dump(prec,info,prefix,head) - implicit none + implicit none type(psb_dprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head - ! len of prefix_ + ! len of prefix_ info = 0 - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = -1 write(psb_err_unit,*) 'Trying to dump a non-built preconditioner' return end if - + call prec%prec%dump(info,prefix,head) - - + + 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 @@ -209,33 +246,33 @@ contains ! 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 + 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 @@ -243,47 +280,47 @@ contains ! 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 + 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 + 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 - implicit none + implicit none type(psb_dprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -303,14 +340,14 @@ contains return 9999 call psb_error_handler(err_act) - + return end subroutine psb_d_precfree subroutine psb_d_prec_free(prec,info) use psb_base_mod - implicit none + implicit none class(psb_dprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -324,7 +361,7 @@ contains me=-1 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%free(info) if (info /= psb_success_) goto 9999 deallocate(prec%prec,stat=info) @@ -343,26 +380,26 @@ contains class(psb_dprec_type), intent(in) :: prec integer(psb_epk_) :: val integer(psb_ipk_) :: i - + val = 0 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then val = val + prec%prec%sizeof() end if - + end function psb_dprec_sizeof subroutine psb_d_prec_clone(prec,precout,info) - implicit none + implicit none class(psb_dprec_type), intent(inout) :: prec class(psb_dprec_type), intent(inout) :: precout integer(psb_ipk_), intent(out) :: info info = psb_success_ call prec%free(info) - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%clone(precout%prec,info) end if - + end subroutine psb_d_prec_clone end module psb_d_prec_type diff --git a/prec/psb_s_bjacprec.f90 b/prec/psb_s_bjacprec.f90 index 2bcf2e02..2cdd184e 100644 --- a/prec/psb_s_bjacprec.f90 +++ b/prec/psb_s_bjacprec.f90 @@ -57,13 +57,11 @@ module psb_s_bjacprec 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 - - character(len=15), parameter, private :: & - & fact_names(0:2)=(/'None ','ILU(n) ',& - & 'ILU(eps) '/) + & fact_names(0:3)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) '/) + private :: psb_s_bjac_sizeof, psb_s_bjac_precdescr, psb_s_bjac_get_nzeros interface subroutine psb_s_bjac_dump(prec,info,prefix,head) diff --git a/prec/psb_s_prec_type.f90 b/prec/psb_s_prec_type.f90 index ded50ca4..f238f407 100644 --- a/prec/psb_s_prec_type.f90 +++ b/prec/psb_s_prec_type.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,8 +27,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Module to define PREC_DATA, !! !! structure for preconditioning. !! @@ -39,7 +39,7 @@ module psb_s_prec_type use psb_s_base_prec_mod type psb_sprec_type - integer(psb_ipk_) :: ictxt + integer(psb_ipk_) :: ictxt class(psb_s_base_prec_type), allocatable :: prec contains procedure, pass(prec) :: psb_s_apply1_vect @@ -54,6 +54,10 @@ 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) :: cseti => psb_scprecseti + procedure, pass(prec) :: csetc => psb_scprecsetc + procedure, pass(prec) :: csetr => psb_scprecsetr + generic, public :: set => cseti, csetc, csetr 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 @@ -102,7 +106,7 @@ module psb_s_prec_type module procedure psb_sprec_sizeof end interface - interface + interface subroutine psb_s_apply2_vect(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_sprec_type, psb_s_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -114,8 +118,8 @@ module psb_s_prec_type real(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_s_apply2_vect end interface - - interface + + interface subroutine psb_s_apply1_vect(prec,x,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_sprec_type, psb_s_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -126,7 +130,7 @@ module psb_s_prec_type real(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_s_apply1_vect end interface - + interface subroutine psb_s_apply2v(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_sprec_type, psb_s_vect_type, psb_spk_ @@ -139,8 +143,8 @@ module psb_s_prec_type real(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_s_apply2v end interface - - interface + + interface subroutine psb_s_apply1v(prec,x,desc_data,info,trans) import :: psb_ipk_, psb_desc_type, psb_sprec_type, psb_s_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -150,56 +154,89 @@ module psb_s_prec_type character(len=1), optional :: trans end subroutine psb_s_apply1v end interface - + + interface + subroutine psb_scprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_sprec_type, psb_sspmat_type, psb_desc_type, psb_spk_, & + & psb_ipk_ + class(psb_sprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_scprecseti + subroutine psb_scprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_sprec_type, psb_sspmat_type, psb_desc_type, psb_spk_, & + & psb_ipk_ + class(psb_sprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_scprecsetr + subroutine psb_scprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + import :: psb_sprec_type, psb_sspmat_type, psb_desc_type, psb_spk_, & + & psb_ipk_ + class(psb_sprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_scprecsetc +end interface + contains subroutine psb_sfile_prec_descr(prec,iout, root) use psb_base_mod - implicit none + implicit none class(psb_sprec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: root integer(psb_ipk_) :: iout_,info - character(len=20) :: name='prec_descr' - - if (present(iout)) then + character(len=20) :: name='prec_descr' + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") end if call prec%prec%descr(iout=iout,root=root) - + end subroutine psb_sfile_prec_descr subroutine psb_s_prec_dump(prec,info,prefix,head) - implicit none + implicit none type(psb_sprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head - ! len of prefix_ + ! len of prefix_ info = 0 - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = -1 write(psb_err_unit,*) 'Trying to dump a non-built preconditioner' return end if - + call prec%prec%dump(info,prefix,head) - - + + 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 @@ -209,33 +246,33 @@ contains ! 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 + 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 @@ -243,47 +280,47 @@ contains ! 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 + 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 + 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 - implicit none + implicit none type(psb_sprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -303,14 +340,14 @@ contains return 9999 call psb_error_handler(err_act) - + return end subroutine psb_s_precfree subroutine psb_s_prec_free(prec,info) use psb_base_mod - implicit none + implicit none class(psb_sprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -324,7 +361,7 @@ contains me=-1 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%free(info) if (info /= psb_success_) goto 9999 deallocate(prec%prec,stat=info) @@ -343,26 +380,26 @@ contains class(psb_sprec_type), intent(in) :: prec integer(psb_epk_) :: val integer(psb_ipk_) :: i - + val = 0 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then val = val + prec%prec%sizeof() end if - + end function psb_sprec_sizeof subroutine psb_s_prec_clone(prec,precout,info) - implicit none + implicit none class(psb_sprec_type), intent(inout) :: prec class(psb_sprec_type), intent(inout) :: precout integer(psb_ipk_), intent(out) :: info info = psb_success_ call prec%free(info) - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%clone(precout%prec,info) end if - + end subroutine psb_s_prec_clone end module psb_s_prec_type diff --git a/prec/psb_z_bjacprec.f90 b/prec/psb_z_bjacprec.f90 index de9b3518..13bb4bcd 100644 --- a/prec/psb_z_bjacprec.f90 +++ b/prec/psb_z_bjacprec.f90 @@ -57,13 +57,11 @@ module psb_z_bjacprec 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 - - character(len=15), parameter, private :: & - & fact_names(0:2)=(/'None ','ILU(n) ',& - & 'ILU(eps) '/) + & fact_names(0:3)=(/'None ','ILU(0) ',& + & 'ILU(n) ','ILU(eps) '/) + private :: psb_z_bjac_sizeof, psb_z_bjac_precdescr, psb_z_bjac_get_nzeros interface subroutine psb_z_bjac_dump(prec,info,prefix,head) diff --git a/prec/psb_z_prec_type.f90 b/prec/psb_z_prec_type.f90 index 8b1a8f4b..4720762c 100644 --- a/prec/psb_z_prec_type.f90 +++ b/prec/psb_z_prec_type.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,8 +27,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Module to define PREC_DATA, !! !! structure for preconditioning. !! @@ -39,7 +39,7 @@ module psb_z_prec_type use psb_z_base_prec_mod type psb_zprec_type - integer(psb_ipk_) :: ictxt + integer(psb_ipk_) :: ictxt class(psb_z_base_prec_type), allocatable :: prec contains procedure, pass(prec) :: psb_z_apply1_vect @@ -54,6 +54,10 @@ 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) :: cseti => psb_zcprecseti + procedure, pass(prec) :: csetc => psb_zcprecsetc + procedure, pass(prec) :: csetr => psb_zcprecsetr + generic, public :: set => cseti, csetc, csetr 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 @@ -102,7 +106,7 @@ module psb_z_prec_type module procedure psb_zprec_sizeof end interface - interface + interface subroutine psb_z_apply2_vect(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_zprec_type, psb_z_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -114,8 +118,8 @@ module psb_z_prec_type complex(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_z_apply2_vect end interface - - interface + + interface subroutine psb_z_apply1_vect(prec,x,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_zprec_type, psb_z_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -126,7 +130,7 @@ module psb_z_prec_type complex(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_z_apply1_vect end interface - + interface subroutine psb_z_apply2v(prec,x,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_zprec_type, psb_z_vect_type, psb_dpk_ @@ -139,8 +143,8 @@ module psb_z_prec_type complex(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_z_apply2v end interface - - interface + + interface subroutine psb_z_apply1v(prec,x,desc_data,info,trans) import :: psb_ipk_, psb_desc_type, psb_zprec_type, psb_z_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -150,56 +154,89 @@ module psb_z_prec_type character(len=1), optional :: trans end subroutine psb_z_apply1v end interface - + + interface + subroutine psb_zcprecseti(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_zprec_type, psb_zspmat_type, psb_desc_type, psb_dpk_, & + & psb_ipk_ + class(psb_zprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_zcprecseti + subroutine psb_zcprecsetr(prec,what,val,info,ilev,ilmax,pos,idx) + import :: psb_zprec_type, psb_zspmat_type, psb_desc_type, psb_dpk_, & + & psb_ipk_ + class(psb_zprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_zcprecsetr + subroutine psb_zcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx) + import :: psb_zprec_type, psb_zspmat_type, psb_desc_type, psb_dpk_, & + & psb_ipk_ + class(psb_zprec_type), intent(inout) :: prec + character(len=*), intent(in) :: what + character(len=*), intent(in) :: string + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev,ilmax,idx + character(len=*), optional, intent(in) :: pos + end subroutine psb_zcprecsetc +end interface + contains subroutine psb_zfile_prec_descr(prec,iout, root) use psb_base_mod - implicit none + implicit none class(psb_zprec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout integer(psb_ipk_), intent(in), optional :: root integer(psb_ipk_) :: iout_,info - character(len=20) :: name='prec_descr' - - if (present(iout)) then + character(len=20) :: name='prec_descr' + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") end if call prec%prec%descr(iout=iout,root=root) - + end subroutine psb_zfile_prec_descr subroutine psb_z_prec_dump(prec,info,prefix,head) - implicit none + implicit none type(psb_zprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head - ! len of prefix_ + ! len of prefix_ info = 0 - if (.not.allocated(prec%prec)) then + if (.not.allocated(prec%prec)) then info = -1 write(psb_err_unit,*) 'Trying to dump a non-built preconditioner' return end if - + call prec%prec%dump(info,prefix,head) - - + + 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 @@ -209,33 +246,33 @@ contains ! 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 + 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 @@ -243,47 +280,47 @@ contains ! 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 + 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 + 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 - implicit none + implicit none type(psb_zprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -303,14 +340,14 @@ contains return 9999 call psb_error_handler(err_act) - + return end subroutine psb_z_precfree subroutine psb_z_prec_free(prec,info) use psb_base_mod - implicit none + implicit none class(psb_zprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: me, err_act,i @@ -324,7 +361,7 @@ contains me=-1 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%free(info) if (info /= psb_success_) goto 9999 deallocate(prec%prec,stat=info) @@ -343,26 +380,26 @@ contains class(psb_zprec_type), intent(in) :: prec integer(psb_epk_) :: val integer(psb_ipk_) :: i - + val = 0 - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then val = val + prec%prec%sizeof() end if - + end function psb_zprec_sizeof subroutine psb_z_prec_clone(prec,precout,info) - implicit none + implicit none class(psb_zprec_type), intent(inout) :: prec class(psb_zprec_type), intent(inout) :: precout integer(psb_ipk_), intent(out) :: info info = psb_success_ call prec%free(info) - if (allocated(prec%prec)) then + if (allocated(prec%prec)) then call prec%prec%clone(precout%prec,info) end if - + end subroutine psb_z_prec_clone end module psb_z_prec_type diff --git a/test/pargen/psb_d_pde2d.f90 b/test/pargen/psb_d_pde2d.f90 index 4b8b6584..9a40bc9d 100644 --- a/test/pargen/psb_d_pde2d.f90 +++ b/test/pargen/psb_d_pde2d.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,23 +27,23 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_d_pde2d.f90 ! ! Program: psb_d_pde2d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 2d ! -! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) +! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) ! - ------ - ------ ----- + ------ + c u = f -! dxdx dydy dx dy +! dxdx dydy dx dy ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit square 0<=x,y<=1. ! @@ -63,31 +63,31 @@ module psb_d_pde2d_mod & psb_dspmat_type, psb_d_vect_type, dzero,& & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type - interface + interface function d_func_2d(x,y) result(val) import :: psb_dpk_ real(psb_dpk_), intent(in) :: x,y real(psb_dpk_) :: val end function d_func_2d - end interface + end interface interface psb_gen_pde2d module procedure psb_d_gen_pde2d end interface psb_gen_pde2d - + contains - + function d_null_func_2d(x,y) result(val) real(psb_dpk_), intent(in) :: x,y real(psb_dpk_) :: val - + val = dzero end function d_null_func_2d ! - ! functions parametrizing the differential equation + ! functions parametrizing the differential equation ! ! @@ -101,48 +101,48 @@ contains ! function b1(x,y) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: b1 real(psb_dpk_), intent(in) :: x,y b1=dzero end function b1 function b2(x,y) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: b2 real(psb_dpk_), intent(in) :: x,y b2=dzero end function b2 function c(x,y) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: c real(psb_dpk_), intent(in) :: x,y c=0.d0 end function c function a1(x,y) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a1 + implicit none + real(psb_dpk_) :: a1 real(psb_dpk_), intent(in) :: x,y a1=done/80 end function a1 function a2(x,y) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: a2 real(psb_dpk_), intent(in) :: x,y a2=done/80 end function a2 function g(x,y) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: g real(psb_dpk_), intent(in) :: x,y g = dzero if (x == done) then g = done - else if (x == dzero) then + else if (x == dzero) then g = exp(-y**2) end if end function g @@ -150,7 +150,7 @@ contains ! ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. + ! the rhs. ! subroutine psb_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,info,& & f,amold,vmold,imold,partition,nrl,iv) @@ -158,13 +158,13 @@ contains use psb_util_mod ! ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) + ! + ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) ! - ------ - ------ + ----- + ------ + c u = f - ! dxdx dydy dx dy + ! dxdx dydy dx dy ! ! with Dirichlet boundary conditions - ! u = g + ! u = g ! ! on the unit square 0<=x,y<=1. ! @@ -221,7 +221,7 @@ contains call psb_info(ictxt, iam, np) - if (present(f)) then + if (present(f)) then f_ => f else f_ => d_null_func_2d @@ -241,9 +241,9 @@ contains else partition_ = 3 end if - + ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes + ! estimate of the number of non zeroes m = (1_psb_lpk_)*idim*idim n = m @@ -252,8 +252,8 @@ contains t0 = psb_wtime() select case(partition_) case(1) - ! A BLOCK partition - if (present(nrl)) then + ! A BLOCK partition + if (present(nrl)) then nr = nrl else ! @@ -264,46 +264,46 @@ contains end if nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) then + call psb_sum(ictxt,nt) + if (nt /= m) then write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! First example of use of CDALL: specify for each process a number of ! contiguous rows - ! + ! call psb_cdall(ictxt,desc_a,info,nl=nr) myidx = desc_a%get_global_indices() nlr = size(myidx) case(2) ! A partition defined by the user through IV - - if (present(iv)) then + + if (present(iv)) then if (size(iv) /= m) then write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if else write(psb_err_unit,*) iam, 'Initialization error: IV not present' info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! + ! process that owns it + ! call psb_cdall(ictxt,desc_a,info,vg=iv) myidx = desc_a%get_global_indices() nlr = size(myidx) @@ -318,7 +318,7 @@ contains npy = npdims(2) allocate(bndx(0:npx),bndy(0:npy)) - ! We can reuse idx2ijk for process indices as well. + ! We can reuse idx2ijk for process indices as well. call idx2ijk(iamx,iamy,iam,npx,npy,base=0) ! Now let's split the 2D square in rectangles call dist1Didx(bndx,idim,npx) @@ -326,7 +326,7 @@ contains call dist1Didx(bndy,idim,npy) myny = bndy(iamy+1)-bndy(iamy) - ! How many indices do I own? + ! How many indices do I own? nlr = mynx*myny allocate(myidx(nlr)) ! Now, let's generate the list of indices I own @@ -348,9 +348,9 @@ contains ! ! Third example of use of CDALL: specify for each process ! the set of global indices it owns. - ! + ! call psb_cdall(ictxt,desc_a,info,vl=myidx) - + case default write(psb_err_unit,*) iam, 'Initialization error: should not get here' info = -1 @@ -359,9 +359,9 @@ contains return end select - + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess + ! define rhs from boundary conditions; also build initial guess if (info == psb_success_) call psb_geall(xv,desc_a,info) if (info == psb_success_) call psb_geall(bv,desc_a,info) @@ -376,12 +376,12 @@ contains end if ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! allocate(val(20*nb),irow(20*nb),& &icol(20*nb),stat=info) - if (info /= psb_success_ ) then + if (info /= psb_success_ ) then info=psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 @@ -394,11 +394,11 @@ contains call psb_barrier(ictxt) t1 = psb_wtime() do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) + ib = min(nb,nlr-ii+1) icoeff = 1 do k=1,ib i=ii+k-1 - ! local matrix pointer + ! local matrix pointer glob_row=myidx(i) ! compute gridpoint coordinates call idx2ijk(ix,iy,glob_row,idim,idim) @@ -408,11 +408,11 @@ contains zt(k) = f_(x,y) ! internal point: build discretization - ! + ! ! term depending on (x-1,y) ! val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2 - if (ix == 1) then + if (ix == 1) then zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix-1,iy,idim,idim) @@ -421,7 +421,7 @@ contains endif ! term depending on (x,y-1) val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2 - if (iy == 1) then + if (iy == 1) then zt(k) = g(x,dzero)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix,iy-1,idim,idim) @@ -433,10 +433,10 @@ contains val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) call ijk2idx(icol(icoeff),ix,iy,idim,idim) irow(icoeff) = glob_row - icoeff = icoeff+1 + icoeff = icoeff+1 ! term depending on (x,y+1) val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2 - if (iy == idim) then + if (iy == idim) then zt(k) = g(x,done)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix,iy+1,idim,idim) @@ -445,7 +445,7 @@ contains endif ! term depending on (x+1,y) val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2 - if (ix==idim) then + if (ix==idim) then zt(k) = g(done,y)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix+1,iy,idim,idim) @@ -479,8 +479,8 @@ contains tcdasb = psb_wtime()-t1 call psb_barrier(ictxt) t1 = psb_wtime() - if (info == psb_success_) then - if (present(amold)) then + if (info == psb_success_) then + if (present(amold)) then call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) else call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) @@ -503,7 +503,7 @@ contains end if tasb = psb_wtime()-t1 call psb_barrier(ictxt) - ttot = psb_wtime() - t0 + ttot = psb_wtime() - t0 call psb_amx(ictxt,talc) call psb_amx(ictxt,tgen) @@ -544,9 +544,9 @@ program psb_d_pde2d integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_dpk_), parameter :: one = done - real(psb_dpk_) :: t1, t2, tprec + real(psb_dpk_) :: t1, t2, tprec ! sparse matrix and preconditioner type(psb_dspmat_type) :: a @@ -563,6 +563,14 @@ program psb_d_pde2d integer(psb_epk_) :: amatsize, precsize, descsize, d2size real(psb_dpk_) :: err, eps + ! Parameters for solvers in Block-Jacobi preconditioner + type ainvparms + character(len=12) :: alg, orth_alg, ilu_alg, ilut_scale + integer(psb_ipk_) :: fill, inv_fill + real(psb_dpk_) :: thresh, inv_thresh + end type ainvparms + type(ainvparms) :: parms + ! other variables integer(psb_ipk_) :: info, i character(len=20) :: name,ch_err @@ -574,7 +582,7 @@ program psb_d_pde2d call psb_init(ictxt) call psb_info(ictxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ictxt) stop @@ -585,21 +593,21 @@ program psb_d_pde2d ! ! Hello world ! - if (iam == psb_root_) then + if (iam == psb_root_) then write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde2d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) + call psb_gen_pde2d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -612,9 +620,28 @@ program psb_d_pde2d if (iam == psb_root_) write(psb_out_unit,'(" ")') ! ! prepare the preconditioner. - ! + ! if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype call prec%init(ictxt,ptype,info) + ! + ! Set the options for the BJAC preconditioner + ! + if (psb_toupper(ptype) == "BJAC") then + call prec%set('sub_solve', parms%alg, info) + select case (psb_toupper(parms%alg)) + case ("ILU") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('ilu_alg', parms%ilu_alg, info) + case ("ILUT") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('sub_iluthrs', parms%thresh, info) + call prec%set('ilut_scale', parms%ilut_scale, info) + case default + ! Do nothing, use default setting in the init routine + end select + else + ! nothing to set for NONE or DIAG preconditioner + end if call psb_barrier(ictxt) t1 = psb_wtime() @@ -634,14 +661,14 @@ program psb_d_pde2d if (iam == psb_root_) write(psb_out_unit,'(" ")') call prec%descr() ! - ! iterative method parameters + ! iterative method parameters ! if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd call psb_barrier(ictxt) - t1 = psb_wtime() + t1 = psb_wtime() eps = 1.d-6 - call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& - & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) + call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& + & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) if(info /= psb_success_) then info=psb_err_from_subroutine_ @@ -671,14 +698,14 @@ program psb_d_pde2d write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err write(psb_out_unit,'("Info on exit : ",i12)')info write(psb_out_unit,'("Total memory occupation for A: ",i12)')amatsize - write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize write(psb_out_unit,'("Storage format for A: ",a)') a%get_fmt() write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt() end if - ! + ! ! cleanup storage and exit ! call psb_gefree(bv,desc_a,info) @@ -704,13 +731,14 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms) integer(psb_ipk_) :: ictxt character(len=*) :: kmethd, ptype, afmt integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst,ipart integer(psb_ipk_) :: np, iam integer(psb_ipk_) :: ip, inp_unit character(len=1024) :: filename + type(ainvparms) :: parms call psb_info(ictxt, iam, np) @@ -739,12 +767,12 @@ contains if (ip >= 4) then read(inp_unit,*) ipart else - ipart = 3 + ipart = 3 endif if (ip >= 5) then read(inp_unit,*) istopc else - istopc=1 + istopc=1 endif if (ip >= 6) then read(inp_unit,*) itmax @@ -761,8 +789,27 @@ contains else irst=1 endif + if (ip >= 9) then + read(inp_unit,*) parms%alg + read(inp_unit,*) parms%ilu_alg + read(inp_unit,*) parms%ilut_scale + read(inp_unit,*) parms%fill + read(inp_unit,*) parms%inv_fill + read(inp_unit,*) parms%thresh + read(inp_unit,*) parms%inv_thresh + read(inp_unit,*) parms%orth_alg + else + parms%alg = 'ILU' ! Block Solver ILU,ILUT,INVK,AINVT,AORTH + parms%ilu_alg = 'NONE' ! If ILU : MILU or NONE othewise ignored + parms%ilut_scale = 'NONE' ! If ILUT: NONE, MAXVAL, DIAG, ARWSUM, ACLSUM, ARCSUM + parms%fill = 0 ! Level of fill for forward factorization + parms%inv_fill = 1 ! Level of fill for inverse factorization (only INVK) + parms%thresh = 1E-1_psb_dpk_ ! Threshold for forward factorization + parms%inv_thresh = 1E-1_psb_dpk_ ! Threshold for inverse factorization + parms%orth_alg = 'LLK' ! What orthogonalization algorithm? + endif - write(psb_out_unit,'("Solving matrix : ell1")') + write(psb_out_unit,'("Solving matrix : ell1")') write(psb_out_unit,'("Grid dimensions : ",i5," x ",i5)')idim,idim write(psb_out_unit,'("Number of processors : ",i0)') np select case(ipart) @@ -775,11 +822,32 @@ contains write(psb_out_unit,'("Unknown data distrbution, defaulting to 2D")') end select write(psb_out_unit,'("Preconditioner : ",a)') ptype + if( psb_toupper(ptype) == "BJAC" ) then + write(psb_out_unit,'("Block subsolver : ",a)') parms%alg + select case (psb_toupper(parms%alg)) + case ('ILU') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("MILU : ",a)') parms%ilu_alg + case ('ILUT') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale + case ('INVK') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh + case ('AINVT','AORTH') + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh + case default + write(psb_out_unit,'("Unknown diagonal solver")') + end select + end if write(psb_out_unit,'("Iterative method : ",a)') kmethd write(psb_out_unit,'(" ")') else ! wrong number of parameter, print an error message and exit - call pr_usage(izero) + call pr_usage(izero) call psb_abort(ictxt) stop 1 endif @@ -803,15 +871,15 @@ contains end subroutine get_parms ! - ! print an error message - ! + ! print an error message + ! subroutine pr_usage(iout) integer(psb_ipk_) :: iout write(iout,*)'incorrect parameter(s) found' write(iout,*)' usage: pde2d90 methd prec dim & - &[ipart istop itmax itrace]' + &[ipart istop itmax itrace]' write(iout,*)' where:' - write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' write(iout,*)' prec : bjac diag none' write(iout,*)' dim number of points along each axis' write(iout,*)' the size of the resulting linear ' @@ -819,11 +887,9 @@ contains write(iout,*)' ipart data partition 1 3 ' write(iout,*)' istop stopping criterion 1, 2 ' write(iout,*)' itmax maximum number of iterations [500] ' - write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' itrace <=0 (no tracing, default) or ' write(iout,*)' >= 1 do tracing every itrace' - write(iout,*)' iterations ' + write(iout,*)' iterations ' end subroutine pr_usage end program psb_d_pde2d - - diff --git a/test/pargen/psb_d_pde3d.f90 b/test/pargen/psb_d_pde3d.f90 index 1a07fffd..4bac87a4 100644 --- a/test/pargen/psb_d_pde3d.f90 +++ b/test/pargen/psb_d_pde3d.f90 @@ -606,7 +606,7 @@ program psb_d_pde3d ! Parameters for solvers in Block-Jacobi preconditioner type ainvparms - character(len=12) :: alg, orth_alg + character(len=12) :: alg, orth_alg, ilu_alg, ilut_scale integer(psb_ipk_) :: fill, inv_fill real(psb_dpk_) :: thresh, inv_thresh end type ainvparms @@ -664,6 +664,25 @@ program psb_d_pde3d ! if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype call prec%init(ictxt,ptype,info) + ! + ! Set the options for the BJAC preconditioner + ! + if (psb_toupper(ptype) == "BJAC") then + call prec%set('sub_solve', parms%alg, info) + select case (psb_toupper(parms%alg)) + case ("ILU") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('ilu_alg', parms%ilu_alg, info) + case ("ILUT") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('sub_iluthrs', parms%thresh, info) + call prec%set('ilut_scale', parms%ilut_scale, info) + case default + ! Do nothing, use default setting in the init routine + end select + else + ! nothing to set for NONE or DIAG preconditioner + end if call psb_barrier(ictxt) t1 = psb_wtime() @@ -813,16 +832,20 @@ contains irst=1 endif if (ip >= 9) then - read(psb_inp_unit,*) parms%alg - read(psb_inp_unit,*) parms%fill - read(psb_inp_unit,*) parms%inv_fill - read(psb_inp_unit,*) parms%thresh - read(psb_inp_unit,*) parms%inv_thresh - read(psb_inp_unit,*) parms%orth_alg + read(inp_unit,*) parms%alg + read(inp_unit,*) parms%ilu_alg + read(inp_unit,*) parms%ilut_scale + read(inp_unit,*) parms%fill + read(inp_unit,*) parms%inv_fill + read(inp_unit,*) parms%thresh + read(inp_unit,*) parms%inv_thresh + read(inp_unit,*) parms%orth_alg else - parms%alg = 'ILU' ! AINV variant: ILU,ILUT,MILU,INVK,AINVT,AORTH - parms%fill = 0 ! Fill in for forward factorization - parms%inv_fill = 1 ! Fill in for inverse factorization + parms%alg = 'ILU' ! Block Solver ILU,ILUT,INVK,AINVT,AORTH + parms%ilu_alg = 'NONE' ! If ILU : MILU or NONE othewise ignored + parms%ilut_scale = 'NONE' ! If ILUT: NONE, MAXVAL, DIAG, ARWSUM, ACLSUM, ARCSUM + parms%fill = 0 ! Level of fill for forward factorization + parms%inv_fill = 1 ! Level of fill for inverse factorization (only INVK) parms%thresh = 1E-1_psb_dpk_ ! Threshold for forward factorization parms%inv_thresh = 1E-1_psb_dpk_ ! Threshold for inverse factorization parms%orth_alg = 'LLK' ! What orthogonalization algorithm? @@ -846,16 +869,20 @@ contains if( psb_toupper(ptype) == "BJAC" ) then write(psb_out_unit,'("Block subsolver : ",a)') parms%alg select case (psb_toupper(parms%alg)) - case ('ILU','ILUT','MILU') + case ('ILU') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("MILU : ",a)') parms%ilu_alg + case ('ILUT') write(psb_out_unit,'("Fill in : ",i0)') parms%fill - write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale case ('INVK') - write(psb_out_unit,'("Fill : ",i0)') parms%fill - write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh - write(psb_out_unit,'("Invese Fill : ",i0)') parms%inv_fill - write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh case ('AINVT','AORTH') - write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh case default write(psb_out_unit,'("Unknown diagonal solver")') end select diff --git a/test/pargen/psb_s_pde2d.f90 b/test/pargen/psb_s_pde2d.f90 index b0fe9a7e..54b04ba7 100644 --- a/test/pargen/psb_s_pde2d.f90 +++ b/test/pargen/psb_s_pde2d.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,23 +27,23 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_s_pde2d.f90 ! ! Program: psb_s_pde2d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 2d ! -! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) +! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) ! - ------ - ------ ----- + ------ + c u = f -! dxdx dydy dx dy +! dxdx dydy dx dy ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit square 0<=x,y<=1. ! @@ -63,31 +63,31 @@ module psb_s_pde2d_mod & psb_sspmat_type, psb_s_vect_type, szero,& & psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type - interface + interface function s_func_2d(x,y) result(val) import :: psb_spk_ real(psb_spk_), intent(in) :: x,y real(psb_spk_) :: val end function s_func_2d - end interface + end interface interface psb_gen_pde2d module procedure psb_s_gen_pde2d end interface psb_gen_pde2d - + contains - + function s_null_func_2d(x,y) result(val) real(psb_spk_), intent(in) :: x,y real(psb_spk_) :: val - + val = szero end function s_null_func_2d ! - ! functions parametrizing the differential equation + ! functions parametrizing the differential equation ! ! @@ -101,48 +101,48 @@ contains ! function b1(x,y) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: b1 real(psb_spk_), intent(in) :: x,y b1=szero end function b1 function b2(x,y) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: b2 real(psb_spk_), intent(in) :: x,y b2=szero end function b2 function c(x,y) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: c real(psb_spk_), intent(in) :: x,y c=0.d0 end function c function a1(x,y) use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: a1 + implicit none + real(psb_spk_) :: a1 real(psb_spk_), intent(in) :: x,y a1=sone/80 end function a1 function a2(x,y) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: a2 real(psb_spk_), intent(in) :: x,y a2=sone/80 end function a2 function g(x,y) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: g real(psb_spk_), intent(in) :: x,y g = szero if (x == sone) then g = sone - else if (x == szero) then + else if (x == szero) then g = exp(-y**2) end if end function g @@ -150,7 +150,7 @@ contains ! ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. + ! the rhs. ! subroutine psb_s_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,info,& & f,amold,vmold,imold,partition,nrl,iv) @@ -158,13 +158,13 @@ contains use psb_util_mod ! ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) + ! + ! a1 dd(u) a2 dd(u) b1 d(u) b2 d(u) ! - ------ - ------ + ----- + ------ + c u = f - ! dxdx dydy dx dy + ! dxdx dydy dx dy ! ! with Dirichlet boundary conditions - ! u = g + ! u = g ! ! on the unit square 0<=x,y<=1. ! @@ -221,7 +221,7 @@ contains call psb_info(ictxt, iam, np) - if (present(f)) then + if (present(f)) then f_ => f else f_ => s_null_func_2d @@ -241,9 +241,9 @@ contains else partition_ = 3 end if - + ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes + ! estimate of the number of non zeroes m = (1_psb_lpk_)*idim*idim n = m @@ -252,8 +252,8 @@ contains t0 = psb_wtime() select case(partition_) case(1) - ! A BLOCK partition - if (present(nrl)) then + ! A BLOCK partition + if (present(nrl)) then nr = nrl else ! @@ -264,46 +264,46 @@ contains end if nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) then + call psb_sum(ictxt,nt) + if (nt /= m) then write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! First example of use of CDALL: specify for each process a number of ! contiguous rows - ! + ! call psb_cdall(ictxt,desc_a,info,nl=nr) myidx = desc_a%get_global_indices() nlr = size(myidx) case(2) ! A partition defined by the user through IV - - if (present(iv)) then + + if (present(iv)) then if (size(iv) /= m) then write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if else write(psb_err_unit,*) iam, 'Initialization error: IV not present' info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! + ! process that owns it + ! call psb_cdall(ictxt,desc_a,info,vg=iv) myidx = desc_a%get_global_indices() nlr = size(myidx) @@ -318,7 +318,7 @@ contains npy = npdims(2) allocate(bndx(0:npx),bndy(0:npy)) - ! We can reuse idx2ijk for process indices as well. + ! We can reuse idx2ijk for process indices as well. call idx2ijk(iamx,iamy,iam,npx,npy,base=0) ! Now let's split the 2D square in rectangles call dist1Didx(bndx,idim,npx) @@ -326,7 +326,7 @@ contains call dist1Didx(bndy,idim,npy) myny = bndy(iamy+1)-bndy(iamy) - ! How many indices do I own? + ! How many indices do I own? nlr = mynx*myny allocate(myidx(nlr)) ! Now, let's generate the list of indices I own @@ -348,9 +348,9 @@ contains ! ! Third example of use of CDALL: specify for each process ! the set of global indices it owns. - ! + ! call psb_cdall(ictxt,desc_a,info,vl=myidx) - + case default write(psb_err_unit,*) iam, 'Initialization error: should not get here' info = -1 @@ -359,9 +359,9 @@ contains return end select - + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess + ! define rhs from boundary conditions; also build initial guess if (info == psb_success_) call psb_geall(xv,desc_a,info) if (info == psb_success_) call psb_geall(bv,desc_a,info) @@ -376,12 +376,12 @@ contains end if ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! allocate(val(20*nb),irow(20*nb),& &icol(20*nb),stat=info) - if (info /= psb_success_ ) then + if (info /= psb_success_ ) then info=psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 @@ -394,11 +394,11 @@ contains call psb_barrier(ictxt) t1 = psb_wtime() do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) + ib = min(nb,nlr-ii+1) icoeff = 1 do k=1,ib i=ii+k-1 - ! local matrix pointer + ! local matrix pointer glob_row=myidx(i) ! compute gridpoint coordinates call idx2ijk(ix,iy,glob_row,idim,idim) @@ -408,11 +408,11 @@ contains zt(k) = f_(x,y) ! internal point: build discretization - ! + ! ! term depending on (x-1,y) ! val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2 - if (ix == 1) then + if (ix == 1) then zt(k) = g(szero,y)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix-1,iy,idim,idim) @@ -421,7 +421,7 @@ contains endif ! term depending on (x,y-1) val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2 - if (iy == 1) then + if (iy == 1) then zt(k) = g(x,szero)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix,iy-1,idim,idim) @@ -433,10 +433,10 @@ contains val(icoeff)=(2*sone)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y) call ijk2idx(icol(icoeff),ix,iy,idim,idim) irow(icoeff) = glob_row - icoeff = icoeff+1 + icoeff = icoeff+1 ! term depending on (x,y+1) val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2 - if (iy == idim) then + if (iy == idim) then zt(k) = g(x,sone)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix,iy+1,idim,idim) @@ -445,7 +445,7 @@ contains endif ! term depending on (x+1,y) val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2 - if (ix==idim) then + if (ix==idim) then zt(k) = g(sone,y)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix+1,iy,idim,idim) @@ -479,8 +479,8 @@ contains tcdasb = psb_wtime()-t1 call psb_barrier(ictxt) t1 = psb_wtime() - if (info == psb_success_) then - if (present(amold)) then + if (info == psb_success_) then + if (present(amold)) then call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) else call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) @@ -503,7 +503,7 @@ contains end if tasb = psb_wtime()-t1 call psb_barrier(ictxt) - ttot = psb_wtime() - t0 + ttot = psb_wtime() - t0 call psb_amx(ictxt,talc) call psb_amx(ictxt,tgen) @@ -544,9 +544,9 @@ program psb_s_pde2d integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_spk_), parameter :: one = sone - real(psb_dpk_) :: t1, t2, tprec + real(psb_dpk_) :: t1, t2, tprec ! sparse matrix and preconditioner type(psb_sspmat_type) :: a @@ -563,6 +563,14 @@ program psb_s_pde2d integer(psb_epk_) :: amatsize, precsize, descsize, d2size real(psb_spk_) :: err, eps + ! Parameters for solvers in Block-Jacobi preconditioner + type ainvparms + character(len=12) :: alg, orth_alg, ilu_alg, ilut_scale + integer(psb_ipk_) :: fill, inv_fill + real(psb_spk_) :: thresh, inv_thresh + end type ainvparms + type(ainvparms) :: parms + ! other variables integer(psb_ipk_) :: info, i character(len=20) :: name,ch_err @@ -574,7 +582,7 @@ program psb_s_pde2d call psb_init(ictxt) call psb_info(ictxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ictxt) stop @@ -585,21 +593,21 @@ program psb_s_pde2d ! ! Hello world ! - if (iam == psb_root_) then + if (iam == psb_root_) then write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde2d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) + call psb_gen_pde2d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -612,9 +620,28 @@ program psb_s_pde2d if (iam == psb_root_) write(psb_out_unit,'(" ")') ! ! prepare the preconditioner. - ! + ! if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype call prec%init(ictxt,ptype,info) + ! + ! Set the options for the BJAC preconditioner + ! + if (psb_toupper(ptype) == "BJAC") then + call prec%set('sub_solve', parms%alg, info) + select case (psb_toupper(parms%alg)) + case ("ILU") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('ilu_alg', parms%ilu_alg, info) + case ("ILUT") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('sub_iluthrs', parms%thresh, info) + call prec%set('ilut_scale', parms%ilut_scale, info) + case default + ! Do nothing, use default setting in the init routine + end select + else + ! nothing to set for NONE or DIAG preconditioner + end if call psb_barrier(ictxt) t1 = psb_wtime() @@ -634,14 +661,14 @@ program psb_s_pde2d if (iam == psb_root_) write(psb_out_unit,'(" ")') call prec%descr() ! - ! iterative method parameters + ! iterative method parameters ! if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd call psb_barrier(ictxt) - t1 = psb_wtime() + t1 = psb_wtime() eps = 1.d-6 - call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& - & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) + call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& + & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) if(info /= psb_success_) then info=psb_err_from_subroutine_ @@ -671,14 +698,14 @@ program psb_s_pde2d write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err write(psb_out_unit,'("Info on exit : ",i12)')info write(psb_out_unit,'("Total memory occupation for A: ",i12)')amatsize - write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize write(psb_out_unit,'("Storage format for A: ",a)') a%get_fmt() write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt() end if - ! + ! ! cleanup storage and exit ! call psb_gefree(bv,desc_a,info) @@ -704,13 +731,14 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms) integer(psb_ipk_) :: ictxt character(len=*) :: kmethd, ptype, afmt integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst,ipart integer(psb_ipk_) :: np, iam integer(psb_ipk_) :: ip, inp_unit character(len=1024) :: filename + type(ainvparms) :: parms call psb_info(ictxt, iam, np) @@ -739,12 +767,12 @@ contains if (ip >= 4) then read(inp_unit,*) ipart else - ipart = 3 + ipart = 3 endif if (ip >= 5) then read(inp_unit,*) istopc else - istopc=1 + istopc=1 endif if (ip >= 6) then read(inp_unit,*) itmax @@ -761,8 +789,27 @@ contains else irst=1 endif + if (ip >= 9) then + read(inp_unit,*) parms%alg + read(inp_unit,*) parms%ilu_alg + read(inp_unit,*) parms%ilut_scale + read(inp_unit,*) parms%fill + read(inp_unit,*) parms%inv_fill + read(inp_unit,*) parms%thresh + read(inp_unit,*) parms%inv_thresh + read(inp_unit,*) parms%orth_alg + else + parms%alg = 'ILU' ! Block Solver ILU,ILUT,INVK,AINVT,AORTH + parms%ilu_alg = 'NONE' ! If ILU : MILU or NONE othewise ignored + parms%ilut_scale = 'NONE' ! If ILUT: NONE, MAXVAL, DIAG, ARWSUM, ACLSUM, ARCSUM + parms%fill = 0 ! Level of fill for forward factorization + parms%inv_fill = 1 ! Level of fill for inverse factorization (only INVK) + parms%thresh = 1E-1_psb_spk_ ! Threshold for forward factorization + parms%inv_thresh = 1E-1_psb_spk_ ! Threshold for inverse factorization + parms%orth_alg = 'LLK' ! What orthogonalization algorithm? + endif - write(psb_out_unit,'("Solving matrix : ell1")') + write(psb_out_unit,'("Solving matrix : ell1")') write(psb_out_unit,'("Grid dimensions : ",i5," x ",i5)')idim,idim write(psb_out_unit,'("Number of processors : ",i0)') np select case(ipart) @@ -775,11 +822,32 @@ contains write(psb_out_unit,'("Unknown data distrbution, defaulting to 2D")') end select write(psb_out_unit,'("Preconditioner : ",a)') ptype + if( psb_toupper(ptype) == "BJAC" ) then + write(psb_out_unit,'("Block subsolver : ",a)') parms%alg + select case (psb_toupper(parms%alg)) + case ('ILU') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("MILU : ",a)') parms%ilu_alg + case ('ILUT') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale + case ('INVK') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh + case ('AINVT','AORTH') + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh + case default + write(psb_out_unit,'("Unknown diagonal solver")') + end select + end if write(psb_out_unit,'("Iterative method : ",a)') kmethd write(psb_out_unit,'(" ")') else ! wrong number of parameter, print an error message and exit - call pr_usage(izero) + call pr_usage(izero) call psb_abort(ictxt) stop 1 endif @@ -803,15 +871,15 @@ contains end subroutine get_parms ! - ! print an error message - ! + ! print an error message + ! subroutine pr_usage(iout) integer(psb_ipk_) :: iout write(iout,*)'incorrect parameter(s) found' write(iout,*)' usage: pde2d90 methd prec dim & - &[ipart istop itmax itrace]' + &[ipart istop itmax itrace]' write(iout,*)' where:' - write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' write(iout,*)' prec : bjac diag none' write(iout,*)' dim number of points along each axis' write(iout,*)' the size of the resulting linear ' @@ -819,11 +887,9 @@ contains write(iout,*)' ipart data partition 1 3 ' write(iout,*)' istop stopping criterion 1, 2 ' write(iout,*)' itmax maximum number of iterations [500] ' - write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' itrace <=0 (no tracing, default) or ' write(iout,*)' >= 1 do tracing every itrace' - write(iout,*)' iterations ' + write(iout,*)' iterations ' end subroutine pr_usage end program psb_s_pde2d - - diff --git a/test/pargen/psb_s_pde3d.f90 b/test/pargen/psb_s_pde3d.f90 index d4ad1492..00f9037d 100644 --- a/test/pargen/psb_s_pde3d.f90 +++ b/test/pargen/psb_s_pde3d.f90 @@ -606,9 +606,9 @@ program psb_s_pde3d ! Parameters for solvers in Block-Jacobi preconditioner type ainvparms - character(len=12) :: alg, orth_alg + character(len=12) :: alg, orth_alg, ilu_alg, ilut_scale integer(psb_ipk_) :: fill, inv_fill - real(psb_dpk_) :: thresh, inv_thresh + real(psb_spk_) :: thresh, inv_thresh end type ainvparms type(ainvparms) :: parms @@ -664,6 +664,25 @@ program psb_s_pde3d ! if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype call prec%init(ictxt,ptype,info) + ! + ! Set the options for the BJAC preconditioner + ! + if (psb_toupper(ptype) == "BJAC") then + call prec%set('sub_solve', parms%alg, info) + select case (psb_toupper(parms%alg)) + case ("ILU") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('ilu_alg', parms%ilu_alg, info) + case ("ILUT") + call prec%set('sub_fillin', parms%fill, info) + call prec%set('sub_iluthrs', parms%thresh, info) + call prec%set('ilut_scale', parms%ilut_scale, info) + case default + ! Do nothing, use default setting in the init routine + end select + else + ! nothing to set for NONE or DIAG preconditioner + end if call psb_barrier(ictxt) t1 = psb_wtime() @@ -813,16 +832,20 @@ contains irst=1 endif if (ip >= 9) then - read(psb_inp_unit,*) parms%alg - read(psb_inp_unit,*) parms%fill - read(psb_inp_unit,*) parms%inv_fill - read(psb_inp_unit,*) parms%thresh - read(psb_inp_unit,*) parms%inv_thresh - read(psb_inp_unit,*) parms%orth_alg + read(inp_unit,*) parms%alg + read(inp_unit,*) parms%ilu_alg + read(inp_unit,*) parms%ilut_scale + read(inp_unit,*) parms%fill + read(inp_unit,*) parms%inv_fill + read(inp_unit,*) parms%thresh + read(inp_unit,*) parms%inv_thresh + read(inp_unit,*) parms%orth_alg else - parms%alg = 'ILU' ! AINV variant: ILU,ILUT,MILU,INVK,AINVT,AORTH - parms%fill = 0 ! Fill in for forward factorization - parms%inv_fill = 1 ! Fill in for inverse factorization + parms%alg = 'ILU' ! Block Solver ILU,ILUT,INVK,AINVT,AORTH + parms%ilu_alg = 'NONE' ! If ILU : MILU or NONE othewise ignored + parms%ilut_scale = 'NONE' ! If ILUT: NONE, MAXVAL, DIAG, ARWSUM, ACLSUM, ARCSUM + parms%fill = 0 ! Level of fill for forward factorization + parms%inv_fill = 1 ! Level of fill for inverse factorization (only INVK) parms%thresh = 1E-1_psb_spk_ ! Threshold for forward factorization parms%inv_thresh = 1E-1_psb_spk_ ! Threshold for inverse factorization parms%orth_alg = 'LLK' ! What orthogonalization algorithm? @@ -846,16 +869,20 @@ contains if( psb_toupper(ptype) == "BJAC" ) then write(psb_out_unit,'("Block subsolver : ",a)') parms%alg select case (psb_toupper(parms%alg)) - case ('ILU','ILUT','MILU') + case ('ILU') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("MILU : ",a)') parms%ilu_alg + case ('ILUT') write(psb_out_unit,'("Fill in : ",i0)') parms%fill - write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale case ('INVK') - write(psb_out_unit,'("Fill : ",i0)') parms%fill - write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh - write(psb_out_unit,'("Invese Fill : ",i0)') parms%inv_fill - write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh + write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh case ('AINVT','AORTH') - write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh case default write(psb_out_unit,'("Unknown diagonal solver")') end select diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index d54ce5b6..f4b45430 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -1,4 +1,4 @@ -8 Number of entries below this +17 Number of entries below this BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES FCG CGR BJAC Preconditioner NONE DIAG BJAC CSR Storage format for matrix A: CSR COO @@ -8,11 +8,11 @@ CSR Storage format for matrix A: CSR COO 0100 MAXIT 05 ITRACE 002 IRST restart for RGMRES and BiCGSTABL -ILU Factorization variant: ILU,ILUT,MILU,INVK,AINVT,AORTH -0 Fill in for forward factorization -1 Fill in for inverse factorization (ignored if not INVK) -1E-1 Threshold for forward factorization (ignored if ILU) -1E-1 Threshold for inverse factorization (ignored if ILU,ILUT,MILU) -LLK What orthogonalization algorithm? (ignored if ILU,ILUT,MILU,INVK) - - +ILU Block Solver ILU,ILUT,INVK,AINVT,AORTH +NONE If ILU : MILU or NONE othewise ignored +NONE Scaling if ILUT: NONE, MAXVAL otherwise ignored +0 Level of fill for forward factorization +1 Level of fill for inverse factorization (only INVK) +1E-1 Threshold for forward factorization +1E-1 Threshold for inverse factorization (Only INVK, AINVT) +LLK What orthogonalization algorithm? (Only AINVT)