diff --git a/mlprec/Makefile b/mlprec/Makefile index 3d3adced..c32d6d4f 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -11,7 +11,7 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBINCDIR) $(FMFLAG)$(PSBLIBD DMODOBJS=mld_d_prec_type.o mld_d_ilu_fact_mod.o \ mld_d_inner_mod.o mld_d_ilu_solver.o mld_d_diag_solver.o mld_d_jac_smoother.o mld_d_as_smoother.o \ mld_d_umf_solver.o mld_d_slu_solver.o mld_d_sludist_solver.o mld_d_id_solver.o\ - mld_d_base_solver_mod.o mld_d_base_smoother_mod.o mld_d_onelev_mod.o + mld_d_base_solver_mod.o mld_d_base_smoother_mod.o mld_d_onelev_mod.o mld_d_gs_solver.o SMODOBJS=mld_s_prec_type.o mld_s_ilu_fact_mod.o \ mld_s_inner_mod.o mld_s_ilu_solver.o mld_s_diag_solver.o mld_s_jac_smoother.o mld_s_as_smoother.o \ @@ -93,7 +93,7 @@ mld_z_base_smoother_mod.o: mld_z_base_solver_mod.o mld_s_base_solver_mod.o mld_d_base_solver_mod.o mld_c_base_solver_mod.o mld_z_base_solver_mod.o: mld_base_prec_type.o -mld_d_id_solver.o mld_d_sludist_solver.o mld_d_slu_solver.o \ +mld_d_gs_solver.o mld_d_id_solver.o mld_d_sludist_solver.o mld_d_slu_solver.o \ mld_d_umf_solver.o mld_d_diag_solver.o mld_d_ilu_solver.o: mld_d_base_solver_mod.o mld_d_prec_type.o mld_d_ilu_fact_mod.o: mld_base_prec_type.o mld_d_base_solver_mod.o diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index 9ffbf8b1..00253aa0 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -161,7 +161,8 @@ module mld_base_prec_type integer(psb_ipk_), parameter :: mld_coarse_subsolve_ = 33 integer(psb_ipk_), parameter :: mld_smoother_sweeps_ = 34 integer(psb_ipk_), parameter :: mld_coarse_aggr_size_ = 35 - integer(psb_ipk_), parameter :: mld_ifpsz_ = 36 + integer(psb_ipk_), parameter :: mld_solver_sweeps_ = 36 + integer(psb_ipk_), parameter :: mld_ifpsz_ = 37 ! ! Legal values for entry: mld_smoother_type_ @@ -275,6 +276,7 @@ module mld_base_prec_type integer(psb_ipk_), parameter :: mld_aggr_thresh_ = 3 integer(psb_ipk_), parameter :: mld_coarse_iluthrs_ = 4 integer(psb_ipk_), parameter :: mld_aggr_scale_ = 5 + integer(psb_ipk_), parameter :: mld_solver_eps_ = 6 integer(psb_ipk_), parameter :: mld_rfpsz_ = 8 ! diff --git a/mlprec/mld_d_gs_solver.f90 b/mlprec/mld_d_gs_solver.f90 new file mode 100644 index 00000000..0d91ecc6 --- /dev/null +++ b/mlprec/mld_d_gs_solver.f90 @@ -0,0 +1,504 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 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 +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ 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 mld_d_gs_solver + + use mld_d_base_solver_mod + + type, extends(mld_d_base_solver_type) :: mld_d_gs_solver_type + type(psb_dspmat_type) :: l, u + integer(psb_ipk_) :: sweeps + real(psb_dpk_) :: eps + contains + procedure, pass(sv) :: dump => mld_d_gs_solver_dmp + procedure, pass(sv) :: ccheck => d_gs_solver_check + procedure, pass(sv) :: clone => mld_d_gs_solver_clone + procedure, pass(sv) :: build => mld_d_gs_solver_bld + procedure, pass(sv) :: cnv => mld_d_gs_solver_cnv + procedure, pass(sv) :: apply_v => mld_d_gs_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_d_gs_solver_apply + procedure, pass(sv) :: free => d_gs_solver_free + procedure, pass(sv) :: seti => d_gs_solver_seti + procedure, pass(sv) :: setc => d_gs_solver_setc + procedure, pass(sv) :: setr => d_gs_solver_setr + procedure, pass(sv) :: cseti => d_gs_solver_cseti + procedure, pass(sv) :: csetc => d_gs_solver_csetc + procedure, pass(sv) :: csetr => d_gs_solver_csetr + procedure, pass(sv) :: descr => d_gs_solver_descr + procedure, pass(sv) :: default => d_gs_solver_default + procedure, pass(sv) :: sizeof => d_gs_solver_sizeof + procedure, pass(sv) :: get_nzeros => d_gs_solver_get_nzeros + procedure, nopass :: get_fmt => d_gs_solver_get_fmt + end type mld_d_gs_solver_type + + + private :: d_gs_solver_bld, d_gs_solver_apply, & + & d_gs_solver_free, d_gs_solver_seti, & + & d_gs_solver_setc, d_gs_solver_setr,& + & d_gs_solver_descr, d_gs_solver_sizeof, & + & d_gs_solver_default, d_gs_solver_dmp, & + & d_gs_solver_apply_vect, d_gs_solver_get_nzeros, & + & d_gs_solver_get_fmt, d_gs_solver_check + + + interface + subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_d_gs_solver_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_gs_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_gs_solver_apply_vect + end interface + + interface + subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_d_gs_solver_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_gs_solver_type), intent(inout) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_gs_solver_apply + end interface + + interface + subroutine mld_d_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + import :: psb_desc_type, mld_d_gs_solver_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_gs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_d_gs_solver_bld + end interface + + interface + subroutine mld_d_gs_solver_cnv(sv,info,amold,vmold,imold) + import :: mld_d_gs_solver_type, psb_dpk_, & + & psb_d_base_sparse_mat, psb_d_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + class(mld_d_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_d_gs_solver_cnv + end interface + + interface + subroutine mld_d_gs_solver_dmp(sv,ictxt,level,info,prefix,head,solver) + import :: psb_desc_type, mld_d_gs_solver_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, & + & psb_ipk_ + implicit none + class(mld_d_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(in) :: ictxt + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver + end subroutine mld_d_gs_solver_dmp + end interface + + interface + subroutine mld_d_gs_solver_clone(sv,svout,info) + import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & + & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, & + & mld_d_base_solver_type, mld_d_gs_solver_type, psb_ipk_ + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + class(mld_d_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_gs_solver_clone + end interface + +contains + + subroutine d_gs_solver_default(sv) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + + sv%sweeps = ione + sv%eps = dzero + + return + end subroutine d_gs_solver_default + + subroutine d_gs_solver_check(sv,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + call mld_check_def(sv%sweeps,& + & 'GS Sweeps',ione,is_int_positive) + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine d_gs_solver_check + + + subroutine d_gs_solver_seti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_seti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(what) + case(mld_solver_sweeps_) + sv%sweeps = val + case default + call sv%mld_d_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_seti + + subroutine d_gs_solver_setc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_setc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setc + + subroutine d_gs_solver_setr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_setr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(what) + case(mld_solver_eps_) + sv%eps = val + case default + call sv%mld_d_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setr + + subroutine d_gs_solver_cseti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + case('SOLVER_SWEEPS') + sv%sweeps = val + case default + call sv%mld_d_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_cseti + + subroutine d_gs_solver_csetc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetc + + subroutine d_gs_solver_csetr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('SOLVER_EPS') + sv%eps = val + case default + call sv%mld_d_base_solver_type%set(what,val,info) + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetr + + subroutine d_gs_solver_free(sv,info) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%l%free() + call sv%u%free() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_free + + subroutine d_gs_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(mld_d_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20), parameter :: name='mld_d_gs_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + endif + + if (sv%eps<=dzero) then + write(iout_,*) ' Gauss-Seidel iterative solver with ',& + & sv%sweeps,' sweeps' + else + write(iout_,*) ' Gauss-Seidel iterative solver with tolerance',& + & sv%eps,' and maxit', sv%sweeps + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_descr + + function d_gs_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(mld_d_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = 0 + val = val + sv%l%get_nzeros() + val = val + sv%u%get_nzeros() + + return + end function d_gs_solver_get_nzeros + + function d_gs_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(mld_d_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = psb_sizeof_int + val = val + sv%l%sizeof() + val = val + sv%u%sizeof() + + return + end function d_gs_solver_sizeof + + function d_gs_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Gauss-Seidel solver" + end function d_gs_solver_get_fmt + +end module mld_d_gs_solver diff --git a/mlprec/mld_d_prec_mod.f90 b/mlprec/mld_d_prec_mod.f90 index 82d4a12e..a9474b96 100644 --- a/mlprec/mld_d_prec_mod.f90 +++ b/mlprec/mld_d_prec_mod.f90 @@ -51,6 +51,7 @@ module mld_d_prec_mod use mld_d_id_solver use mld_d_diag_solver use mld_d_ilu_solver + use mld_d_gs_solver interface mld_precinit subroutine mld_dprecinit(p,ptype,info,nlev)