From a2a1532b908feadf2460e1c7e70075a1738c75a8 Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Mon, 30 Nov 2020 23:07:52 +0100 Subject: [PATCH] Added rkr solver --- amgprec/Makefile | 17 +- amgprec/amg_c_rkr_solver.f90 | 590 ++++++++++++++++++ amgprec/amg_d_rkr_solver.f90 | 590 ++++++++++++++++++ amgprec/amg_s_rkr_solver.f90 | 590 ++++++++++++++++++ amgprec/amg_z_rkr_solver.f90 | 590 ++++++++++++++++++ amgprec/impl/solver/Makefile | 6 +- amgprec/impl/solver/amg_c_rkr_solver_impl.f90 | 271 ++++++++ amgprec/impl/solver/amg_d_rkr_solver_impl.f90 | 271 ++++++++ amgprec/impl/solver/amg_s_rkr_solver_impl.f90 | 271 ++++++++ amgprec/impl/solver/amg_z_rkr_solver_impl.f90 | 271 ++++++++ tests/pdegen/amg_d_pde2d.f90 | 117 +++- tests/pdegen/amg_d_pde3d.f90 | 118 +++- tests/pdegen/amg_s_pde2d.f90 | 117 +++- tests/pdegen/amg_s_pde3d.f90 | 118 +++- tests/pdegen/runs/amg_pde2d.inp | 23 +- tests/pdegen/runs/amg_pde3d.inp | 23 +- 16 files changed, 3930 insertions(+), 53 deletions(-) create mode 100644 amgprec/amg_c_rkr_solver.f90 create mode 100644 amgprec/amg_d_rkr_solver.f90 create mode 100644 amgprec/amg_s_rkr_solver.f90 create mode 100644 amgprec/amg_z_rkr_solver.f90 create mode 100644 amgprec/impl/solver/amg_c_rkr_solver_impl.f90 create mode 100644 amgprec/impl/solver/amg_d_rkr_solver_impl.f90 create mode 100644 amgprec/impl/solver/amg_s_rkr_solver_impl.f90 create mode 100644 amgprec/impl/solver/amg_z_rkr_solver_impl.f90 diff --git a/amgprec/Makefile b/amgprec/Makefile index d6f3e631..1c50887a 100644 --- a/amgprec/Makefile +++ b/amgprec/Makefile @@ -15,7 +15,8 @@ DMODOBJS=amg_d_prec_type.o \ amg_d_base_aggregator_mod.o \ amg_d_dec_aggregator_mod.o amg_d_symdec_aggregator_mod.o \ amg_d_ainv_solver.o amg_d_base_ainv_mod.o \ - amg_d_invk_solver.o amg_d_invt_solver.o + amg_d_invk_solver.o amg_d_invt_solver.o \ + amg_d_rkr_solver.o #amg_d_bcmatch_aggregator_mod.o SMODOBJS=amg_s_prec_type.o amg_s_ilu_fact_mod.o \ @@ -26,7 +27,8 @@ SMODOBJS=amg_s_prec_type.o amg_s_ilu_fact_mod.o \ amg_s_base_aggregator_mod.o \ amg_s_dec_aggregator_mod.o amg_s_symdec_aggregator_mod.o \ amg_s_ainv_solver.o amg_s_base_ainv_mod.o \ - amg_s_invk_solver.o amg_s_invt_solver.o + amg_s_invk_solver.o amg_s_invt_solver.o \ + amg_s_rkr_solver.o ZMODOBJS=amg_z_prec_type.o amg_z_ilu_fact_mod.o \ amg_z_inner_mod.o amg_z_ilu_solver.o amg_z_diag_solver.o amg_z_jac_smoother.o amg_z_as_smoother.o \ @@ -36,7 +38,8 @@ ZMODOBJS=amg_z_prec_type.o amg_z_ilu_fact_mod.o \ amg_z_base_aggregator_mod.o \ amg_z_dec_aggregator_mod.o amg_z_symdec_aggregator_mod.o \ amg_z_ainv_solver.o amg_z_base_ainv_mod.o \ - amg_z_invk_solver.o amg_z_invt_solver.o + amg_z_invk_solver.o amg_z_invt_solver.o \ + amg_z_rkr_solver.o CMODOBJS=amg_c_prec_type.o amg_c_ilu_fact_mod.o \ amg_c_inner_mod.o amg_c_ilu_solver.o amg_c_diag_solver.o amg_c_jac_smoother.o amg_c_as_smoother.o \ @@ -46,7 +49,8 @@ CMODOBJS=amg_c_prec_type.o amg_c_ilu_fact_mod.o \ amg_c_base_aggregator_mod.o \ amg_c_dec_aggregator_mod.o amg_c_symdec_aggregator_mod.o \ amg_c_ainv_solver.o amg_c_base_ainv_mod.o \ - amg_c_invk_solver.o amg_c_invt_solver.o + amg_c_invk_solver.o amg_c_invt_solver.o \ + amg_c_rkr_solver.o @@ -137,6 +141,11 @@ amg_c_base_ainv_mod.o: amg_c_base_solver_mod.o amg_base_ainv_mod.o amg_d_base_ainv_mod.o: amg_d_base_solver_mod.o amg_base_ainv_mod.o amg_z_base_ainv_mod.o: amg_z_base_solver_mod.o amg_base_ainv_mod.o +amg_d_rkr_solver.o: amg_d_base_solver_mod.o amg_d_prec_type.o +amg_s_rkr_solver.o: amg_s_base_solver_mod.o amg_s_prec_type.o +amg_c_rkr_solver.o: amg_c_base_solver_mod.o amg_c_prec_type.o +amg_z_rkr_solver.o: amg_z_base_solver_mod.o amg_z_prec_type.o + amg_s_base_solver_mod.o amg_d_base_solver_mod.o amg_c_base_solver_mod.o amg_z_base_solver_mod.o: amg_base_prec_type.o amg_d_mumps_solver.o amg_d_gs_solver.o amg_d_id_solver.o amg_d_sludist_solver.o amg_d_slu_solver.o \ diff --git a/amgprec/amg_c_rkr_solver.f90 b/amgprec/amg_c_rkr_solver.f90 new file mode 100644 index 00000000..fff0756e --- /dev/null +++ b/amgprec/amg_c_rkr_solver.f90 @@ -0,0 +1,590 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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. +! +! +! +! File: amg_c_rkr_solver_mod.f90 +! +! Module: amg_c_rkr_solver_mod +! +module amg_c_rkr_solver + + use amg_c_base_solver_mod + use amg_c_prec_type + + type, extends(amg_c_base_solver_type) :: amg_c_rkr_solver_type + ! + logical :: global + character(len=16) :: method, kprec, sub_solve + integer(psb_ipk_) :: irst, istopc, itmax, itrace, i_sub_solve + integer(psb_ipk_) :: fillin + real(psb_spk_) :: eps + type(amg_cprec_type) :: prec + type(psb_desc_type) :: desc_local + type(psb_cspmat_type) :: a_local + type(psb_c_vect_type) :: x_local, z_local + type(psb_cspmat_type), pointer :: a=>null() + contains + ! + ! + procedure, pass(sv) :: dump => c_rkr_solver_dmp + procedure, pass(sv) :: check => c_rkr_solver_check + procedure, pass(sv) :: clone => c_rkr_solver_clone + procedure, pass(sv) :: clone_settings => c_rkr_solver_clone_settings + procedure, pass(sv) :: cnv => c_rkr_solver_cnv + procedure, pass(sv) :: apply_v => amg_c_rkr_solver_apply_vect + procedure, pass(sv) :: apply_a => amg_c_rkr_solver_apply + procedure, pass(sv) :: clear_data => c_rkr_solver_clear_data + procedure, pass(sv) :: free => c_rkr_solver_free + procedure, pass(sv) :: cseti => c_rkr_solver_cseti + procedure, pass(sv) :: csetc => c_rkr_solver_csetc + procedure, pass(sv) :: csetr => c_rkr_solver_csetr + procedure, pass(sv) :: sizeof => c_rkr_solver_sizeof + procedure, pass(sv) :: get_nzeros => c_rkr_solver_get_nzeros + !procedure, nopass :: get_id => c_rkr_solver_get_id + procedure, pass(sv) :: is_global => c_rkr_solver_is_global + procedure, nopass :: is_iterative => c_rkr_solver_is_iterative + + + ! + ! These methods are specific for the new solver type + ! and therefore need to be overridden + ! + procedure, pass(sv) :: descr => c_rkr_solver_descr + procedure, pass(sv) :: default => c_rkr_solver_default + procedure, pass(sv) :: build => amg_c_rkr_solver_bld + procedure, nopass :: get_fmt => c_rkr_solver_get_fmt + end type amg_c_rkr_solver_type + + + private :: c_rkr_solver_get_fmt, c_rkr_solver_descr, c_rkr_solver_default + + interface + subroutine amg_c_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + import :: psb_desc_type, amg_c_rkr_solver_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_rkr_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_c_vect_type),intent(inout), optional :: initu + end subroutine amg_c_rkr_solver_apply_vect + end interface + + interface + subroutine amg_c_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + import :: psb_desc_type, amg_c_rkr_solver_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_rkr_solver_type), intent(inout) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_spk_),intent(inout), optional :: initu(:) + end subroutine amg_c_rkr_solver_apply + end interface + + interface + subroutine amg_c_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + import :: psb_desc_type, amg_c_rkr_solver_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_c_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_c_rkr_solver_bld + end interface + + +contains + + ! + ! + subroutine c_rkr_solver_default(sv) + + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + + sv%method = 'bicgstab' + sv%kprec = 'bjac' + sv%sub_solve = 'ilu' + sv%i_sub_solve = -1 + sv%fillin = 1 + sv%irst = 30 + sv%istopc = 2 + sv%itmax = 40 + sv%itrace = -1 + sv%eps = 1.d-6 + sv%global = .false. + + return + end subroutine c_rkr_solver_default + + function c_rkr_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(amg_c_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%get_nzeros() + + return + end function c_rkr_solver_get_nzeros + + function c_rkr_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(amg_c_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%sizeof() + sv%desc_local%sizeof() + sv%a_local%sizeof() + + return + end function c_rkr_solver_sizeof + + + subroutine c_rkr_solver_check(sv,info) + + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='c_rkr_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine c_rkr_solver_check + + subroutine c_rkr_solver_cseti(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_c_rkr_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_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='c_rkr_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(trim(what))) + case('RKR_IRST') + sv%irst = val + case('RKR_ISTOPC') + sv%istopc = val + case('RKR_ITMAX') + sv%itmax = val + case('RKR_ITRACE') + sv%itrace = val + case('RKR_SUB_SOLVE') + sv%i_sub_solve = val + case('RKR_FILLIN') + sv%fillin = val + case default + call sv%amg_c_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine c_rkr_solver_cseti + + subroutine c_rkr_solver_csetc(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='c_rkr_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + select case(psb_toupper(trim(what))) + case('RKR_METHOD') + sv%method = psb_toupper(trim(val)) + case('RKR_KPREC') + sv%kprec = psb_toupper(trim(val)) + case('RKR_SUB_SOLVE') + sv%sub_solve = psb_toupper(trim(val)) + case('RKR_GLOBAL') + select case(psb_toupper(trim(val))) + case('LOCAL','FALSE') + sv%global = .false. + case('GLOBAL','TRUE') + sv%global =.true. + end select + case default + call sv%amg_c_base_solver_type%set(what,val,info,idx=idx) + end select + + + 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 c_rkr_solver_csetc + + subroutine c_rkr_solver_csetr(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='c_rkr_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('RKR_EPS') + sv%eps = val + case default + call sv%amg_c_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine c_rkr_solver_csetr + + subroutine c_rkr_solver_clear_data(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='c_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%prec%free(info) + call sv%a_local%free() + call sv%x_local%free(info) + call sv%z_local%free(info) + if ((.not.sv%global).and.sv%desc_local%is_ok()) then + l_ctxt=sv%desc_local%get_context() + call psb_exit(l_ctxt,close=.false.) + end if + call sv%desc_local%free(info) + nullify(sv%a) + call psb_erractionrestore(err_act) + return + end subroutine c_rkr_solver_clear_data + + + subroutine c_rkr_solver_free(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='c_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%clear_data(info) + + call psb_erractionrestore(err_act) + return + end subroutine c_rkr_solver_free + + function c_rkr_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "RKR solver" + end function c_rkr_solver_get_fmt + + subroutine c_rkr_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(amg_c_rkr_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='amg_c_rkr_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = psb_out_unit + endif + + if (sv%global) then + write(iout_,*) ' Recursive Krylov solver (global)' + else + write(iout_,*) ' Recursive Krylov solver (local) ' + end if + write(iout_,*) ' method: ',sv%method + write(iout_,*) ' kprec: ',sv%kprec + if (sv%i_sub_solve > 0) then + write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve) + else + write(iout_,*) ' sub_solve: ',sv%sub_solve + end if + write(iout_,*) ' itmax: ',sv%itmax + write(iout_,*) ' eps: ',sv%eps + write(iout_,*) ' fillin: ',sv%fillin + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine c_rkr_solver_descr + + subroutine c_rkr_solver_cnv(sv,info,amold,vmold,imold) + implicit none + class(amg_c_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + + call sv%prec%cnv(info,amold=amold,vmold=vmold,imold=imold) + + end subroutine c_rkr_solver_cnv + + subroutine c_rkr_solver_clone(sv,svout,info) + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + class(amg_c_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + call svout%free(info) + allocate(svout,stat=info,mold=sv) + select type(so=>svout) + class is(amg_c_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + call sv%a_local%clone(so%a_local,info) + call sv%desc_local%clone(so%desc_local,info) + call sv%prec%clone(so%prec,info) + class default + info = psb_err_internal_error_ + end select + + end subroutine c_rkr_solver_clone + + + subroutine c_rkr_solver_clone_settings(sv,svout,info) + Implicit None + + ! Arguments + class(amg_c_rkr_solver_type), intent(inout) :: sv + class(amg_c_base_solver_type), intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + select type(so=>svout) + class is(amg_c_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + class default + info = psb_err_internal_error_ + end select + + end subroutine c_rkr_solver_clone_settings + + subroutine c_rkr_solver_dmp(sv,desc,level,info,prefix,head,solver,global_num) + implicit none + class(amg_c_rkr_solver_type), intent(in) :: sv + type(psb_desc_type), intent(in) :: desc + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver, global_num + + + call sv%prec%dump(info,prefix=prefix,head=head) + + end subroutine c_rkr_solver_dmp + ! + ! Notify whether RKR is used as a global solver + ! + function c_rkr_solver_is_global(sv) result(val) + implicit none + class(amg_c_rkr_solver_type), intent(in) :: sv + logical :: val + + val = (sv%global) + end function c_rkr_solver_is_global + ! + function c_rkr_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .true. + end function c_rkr_solver_is_iterative + +end module amg_c_rkr_solver diff --git a/amgprec/amg_d_rkr_solver.f90 b/amgprec/amg_d_rkr_solver.f90 new file mode 100644 index 00000000..8405ae72 --- /dev/null +++ b/amgprec/amg_d_rkr_solver.f90 @@ -0,0 +1,590 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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. +! +! +! +! File: amg_d_rkr_solver_mod.f90 +! +! Module: amg_d_rkr_solver_mod +! +module amg_d_rkr_solver + + use amg_d_base_solver_mod + use amg_d_prec_type + + type, extends(amg_d_base_solver_type) :: amg_d_rkr_solver_type + ! + logical :: global + character(len=16) :: method, kprec, sub_solve + integer(psb_ipk_) :: irst, istopc, itmax, itrace, i_sub_solve + integer(psb_ipk_) :: fillin + real(psb_dpk_) :: eps + type(amg_dprec_type) :: prec + type(psb_desc_type) :: desc_local + type(psb_dspmat_type) :: a_local + type(psb_d_vect_type) :: x_local, z_local + type(psb_dspmat_type), pointer :: a=>null() + contains + ! + ! + procedure, pass(sv) :: dump => d_rkr_solver_dmp + procedure, pass(sv) :: check => d_rkr_solver_check + procedure, pass(sv) :: clone => d_rkr_solver_clone + procedure, pass(sv) :: clone_settings => d_rkr_solver_clone_settings + procedure, pass(sv) :: cnv => d_rkr_solver_cnv + procedure, pass(sv) :: apply_v => amg_d_rkr_solver_apply_vect + procedure, pass(sv) :: apply_a => amg_d_rkr_solver_apply + procedure, pass(sv) :: clear_data => d_rkr_solver_clear_data + procedure, pass(sv) :: free => d_rkr_solver_free + procedure, pass(sv) :: cseti => d_rkr_solver_cseti + procedure, pass(sv) :: csetc => d_rkr_solver_csetc + procedure, pass(sv) :: csetr => d_rkr_solver_csetr + procedure, pass(sv) :: sizeof => d_rkr_solver_sizeof + procedure, pass(sv) :: get_nzeros => d_rkr_solver_get_nzeros + !procedure, nopass :: get_id => d_rkr_solver_get_id + procedure, pass(sv) :: is_global => d_rkr_solver_is_global + procedure, nopass :: is_iterative => d_rkr_solver_is_iterative + + + ! + ! These methods are specific for the new solver type + ! and therefore need to be overridden + ! + procedure, pass(sv) :: descr => d_rkr_solver_descr + procedure, pass(sv) :: default => d_rkr_solver_default + procedure, pass(sv) :: build => amg_d_rkr_solver_bld + procedure, nopass :: get_fmt => d_rkr_solver_get_fmt + end type amg_d_rkr_solver_type + + + private :: d_rkr_solver_get_fmt, d_rkr_solver_descr, d_rkr_solver_default + + interface + subroutine amg_d_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + import :: psb_desc_type, amg_d_rkr_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(amg_d_rkr_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(:) + type(psb_d_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_d_vect_type),intent(inout), optional :: initu + end subroutine amg_d_rkr_solver_apply_vect + end interface + + interface + subroutine amg_d_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + import :: psb_desc_type, amg_d_rkr_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(amg_d_rkr_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 + character, intent(in), optional :: init + real(psb_dpk_),intent(inout), optional :: initu(:) + end subroutine amg_d_rkr_solver_apply + end interface + + interface + subroutine amg_d_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + import :: psb_desc_type, amg_d_rkr_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(inout) :: desc_a + class(amg_d_rkr_solver_type), intent(inout) :: sv + 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 amg_d_rkr_solver_bld + end interface + + +contains + + ! + ! + subroutine d_rkr_solver_default(sv) + + Implicit None + + ! Arguments + class(amg_d_rkr_solver_type), intent(inout) :: sv + + sv%method = 'bicgstab' + sv%kprec = 'bjac' + sv%sub_solve = 'ilu' + sv%i_sub_solve = -1 + sv%fillin = 1 + sv%irst = 30 + sv%istopc = 2 + sv%itmax = 40 + sv%itrace = -1 + sv%eps = 1.d-6 + sv%global = .false. + + return + end subroutine d_rkr_solver_default + + function d_rkr_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(amg_d_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%get_nzeros() + + return + end function d_rkr_solver_get_nzeros + + function d_rkr_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(amg_d_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%sizeof() + sv%desc_local%sizeof() + sv%a_local%sizeof() + + return + end function d_rkr_solver_sizeof + + + subroutine d_rkr_solver_check(sv,info) + + Implicit None + + ! Arguments + class(amg_d_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_rkr_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine d_rkr_solver_check + + subroutine d_rkr_solver_cseti(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_d_rkr_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_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_rkr_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(trim(what))) + case('RKR_IRST') + sv%irst = val + case('RKR_ISTOPC') + sv%istopc = val + case('RKR_ITMAX') + sv%itmax = val + case('RKR_ITRACE') + sv%itrace = val + case('RKR_SUB_SOLVE') + sv%i_sub_solve = val + case('RKR_FILLIN') + sv%fillin = val + case default + call sv%amg_d_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_rkr_solver_cseti + + subroutine d_rkr_solver_csetc(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_d_rkr_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_rkr_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + select case(psb_toupper(trim(what))) + case('RKR_METHOD') + sv%method = psb_toupper(trim(val)) + case('RKR_KPREC') + sv%kprec = psb_toupper(trim(val)) + case('RKR_SUB_SOLVE') + sv%sub_solve = psb_toupper(trim(val)) + case('RKR_GLOBAL') + select case(psb_toupper(trim(val))) + case('LOCAL','FALSE') + sv%global = .false. + case('GLOBAL','TRUE') + sv%global =.true. + end select + case default + call sv%amg_d_base_solver_type%set(what,val,info,idx=idx) + end select + + + 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_rkr_solver_csetc + + subroutine d_rkr_solver_csetr(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_d_rkr_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_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_rkr_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('RKR_EPS') + sv%eps = val + case default + call sv%amg_d_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_rkr_solver_csetr + + subroutine d_rkr_solver_clear_data(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_d_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='d_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%prec%free(info) + call sv%a_local%free() + call sv%x_local%free(info) + call sv%z_local%free(info) + if ((.not.sv%global).and.sv%desc_local%is_ok()) then + l_ctxt=sv%desc_local%get_context() + call psb_exit(l_ctxt,close=.false.) + end if + call sv%desc_local%free(info) + nullify(sv%a) + call psb_erractionrestore(err_act) + return + end subroutine d_rkr_solver_clear_data + + + subroutine d_rkr_solver_free(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_d_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='d_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%clear_data(info) + + call psb_erractionrestore(err_act) + return + end subroutine d_rkr_solver_free + + function d_rkr_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "RKR solver" + end function d_rkr_solver_get_fmt + + subroutine d_rkr_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(amg_d_rkr_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='amg_d_rkr_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = psb_out_unit + endif + + if (sv%global) then + write(iout_,*) ' Recursive Krylov solver (global)' + else + write(iout_,*) ' Recursive Krylov solver (local) ' + end if + write(iout_,*) ' method: ',sv%method + write(iout_,*) ' kprec: ',sv%kprec + if (sv%i_sub_solve > 0) then + write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve) + else + write(iout_,*) ' sub_solve: ',sv%sub_solve + end if + write(iout_,*) ' itmax: ',sv%itmax + write(iout_,*) ' eps: ',sv%eps + write(iout_,*) ' fillin: ',sv%fillin + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_rkr_solver_descr + + subroutine d_rkr_solver_cnv(sv,info,amold,vmold,imold) + implicit none + class(amg_d_rkr_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 + + call sv%prec%cnv(info,amold=amold,vmold=vmold,imold=imold) + + end subroutine d_rkr_solver_cnv + + subroutine d_rkr_solver_clone(sv,svout,info) + Implicit None + + ! Arguments + class(amg_d_rkr_solver_type), intent(inout) :: sv + class(amg_d_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + call svout%free(info) + allocate(svout,stat=info,mold=sv) + select type(so=>svout) + class is(amg_d_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + call sv%a_local%clone(so%a_local,info) + call sv%desc_local%clone(so%desc_local,info) + call sv%prec%clone(so%prec,info) + class default + info = psb_err_internal_error_ + end select + + end subroutine d_rkr_solver_clone + + + subroutine d_rkr_solver_clone_settings(sv,svout,info) + Implicit None + + ! Arguments + class(amg_d_rkr_solver_type), intent(inout) :: sv + class(amg_d_base_solver_type), intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + select type(so=>svout) + class is(amg_d_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + class default + info = psb_err_internal_error_ + end select + + end subroutine d_rkr_solver_clone_settings + + subroutine d_rkr_solver_dmp(sv,desc,level,info,prefix,head,solver,global_num) + implicit none + class(amg_d_rkr_solver_type), intent(in) :: sv + type(psb_desc_type), intent(in) :: desc + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver, global_num + + + call sv%prec%dump(info,prefix=prefix,head=head) + + end subroutine d_rkr_solver_dmp + ! + ! Notify whether RKR is used as a global solver + ! + function d_rkr_solver_is_global(sv) result(val) + implicit none + class(amg_d_rkr_solver_type), intent(in) :: sv + logical :: val + + val = (sv%global) + end function d_rkr_solver_is_global + ! + function d_rkr_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .true. + end function d_rkr_solver_is_iterative + +end module amg_d_rkr_solver diff --git a/amgprec/amg_s_rkr_solver.f90 b/amgprec/amg_s_rkr_solver.f90 new file mode 100644 index 00000000..efe60cd2 --- /dev/null +++ b/amgprec/amg_s_rkr_solver.f90 @@ -0,0 +1,590 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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. +! +! +! +! File: amg_s_rkr_solver_mod.f90 +! +! Module: amg_s_rkr_solver_mod +! +module amg_s_rkr_solver + + use amg_s_base_solver_mod + use amg_s_prec_type + + type, extends(amg_s_base_solver_type) :: amg_s_rkr_solver_type + ! + logical :: global + character(len=16) :: method, kprec, sub_solve + integer(psb_ipk_) :: irst, istopc, itmax, itrace, i_sub_solve + integer(psb_ipk_) :: fillin + real(psb_spk_) :: eps + type(amg_sprec_type) :: prec + type(psb_desc_type) :: desc_local + type(psb_sspmat_type) :: a_local + type(psb_s_vect_type) :: x_local, z_local + type(psb_sspmat_type), pointer :: a=>null() + contains + ! + ! + procedure, pass(sv) :: dump => s_rkr_solver_dmp + procedure, pass(sv) :: check => s_rkr_solver_check + procedure, pass(sv) :: clone => s_rkr_solver_clone + procedure, pass(sv) :: clone_settings => s_rkr_solver_clone_settings + procedure, pass(sv) :: cnv => s_rkr_solver_cnv + procedure, pass(sv) :: apply_v => amg_s_rkr_solver_apply_vect + procedure, pass(sv) :: apply_a => amg_s_rkr_solver_apply + procedure, pass(sv) :: clear_data => s_rkr_solver_clear_data + procedure, pass(sv) :: free => s_rkr_solver_free + procedure, pass(sv) :: cseti => s_rkr_solver_cseti + procedure, pass(sv) :: csetc => s_rkr_solver_csetc + procedure, pass(sv) :: csetr => s_rkr_solver_csetr + procedure, pass(sv) :: sizeof => s_rkr_solver_sizeof + procedure, pass(sv) :: get_nzeros => s_rkr_solver_get_nzeros + !procedure, nopass :: get_id => s_rkr_solver_get_id + procedure, pass(sv) :: is_global => s_rkr_solver_is_global + procedure, nopass :: is_iterative => s_rkr_solver_is_iterative + + + ! + ! These methods are specific for the new solver type + ! and therefore need to be overridden + ! + procedure, pass(sv) :: descr => s_rkr_solver_descr + procedure, pass(sv) :: default => s_rkr_solver_default + procedure, pass(sv) :: build => amg_s_rkr_solver_bld + procedure, nopass :: get_fmt => s_rkr_solver_get_fmt + end type amg_s_rkr_solver_type + + + private :: s_rkr_solver_get_fmt, s_rkr_solver_descr, s_rkr_solver_default + + interface + subroutine amg_s_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + import :: psb_desc_type, amg_s_rkr_solver_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_rkr_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_s_vect_type),intent(inout), optional :: initu + end subroutine amg_s_rkr_solver_apply_vect + end interface + + interface + subroutine amg_s_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + import :: psb_desc_type, amg_s_rkr_solver_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_rkr_solver_type), intent(inout) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_spk_),intent(inout), optional :: initu(:) + end subroutine amg_s_rkr_solver_apply + end interface + + interface + subroutine amg_s_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + import :: psb_desc_type, amg_s_rkr_solver_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_s_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_s_rkr_solver_bld + end interface + + +contains + + ! + ! + subroutine s_rkr_solver_default(sv) + + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + + sv%method = 'bicgstab' + sv%kprec = 'bjac' + sv%sub_solve = 'ilu' + sv%i_sub_solve = -1 + sv%fillin = 1 + sv%irst = 30 + sv%istopc = 2 + sv%itmax = 40 + sv%itrace = -1 + sv%eps = 1.d-6 + sv%global = .false. + + return + end subroutine s_rkr_solver_default + + function s_rkr_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(amg_s_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%get_nzeros() + + return + end function s_rkr_solver_get_nzeros + + function s_rkr_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(amg_s_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%sizeof() + sv%desc_local%sizeof() + sv%a_local%sizeof() + + return + end function s_rkr_solver_sizeof + + + subroutine s_rkr_solver_check(sv,info) + + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='s_rkr_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine s_rkr_solver_check + + subroutine s_rkr_solver_cseti(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_s_rkr_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_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='s_rkr_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(trim(what))) + case('RKR_IRST') + sv%irst = val + case('RKR_ISTOPC') + sv%istopc = val + case('RKR_ITMAX') + sv%itmax = val + case('RKR_ITRACE') + sv%itrace = val + case('RKR_SUB_SOLVE') + sv%i_sub_solve = val + case('RKR_FILLIN') + sv%fillin = val + case default + call sv%amg_s_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine s_rkr_solver_cseti + + subroutine s_rkr_solver_csetc(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='s_rkr_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + select case(psb_toupper(trim(what))) + case('RKR_METHOD') + sv%method = psb_toupper(trim(val)) + case('RKR_KPREC') + sv%kprec = psb_toupper(trim(val)) + case('RKR_SUB_SOLVE') + sv%sub_solve = psb_toupper(trim(val)) + case('RKR_GLOBAL') + select case(psb_toupper(trim(val))) + case('LOCAL','FALSE') + sv%global = .false. + case('GLOBAL','TRUE') + sv%global =.true. + end select + case default + call sv%amg_s_base_solver_type%set(what,val,info,idx=idx) + end select + + + 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 s_rkr_solver_csetc + + subroutine s_rkr_solver_csetr(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='s_rkr_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('RKR_EPS') + sv%eps = val + case default + call sv%amg_s_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine s_rkr_solver_csetr + + subroutine s_rkr_solver_clear_data(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='s_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%prec%free(info) + call sv%a_local%free() + call sv%x_local%free(info) + call sv%z_local%free(info) + if ((.not.sv%global).and.sv%desc_local%is_ok()) then + l_ctxt=sv%desc_local%get_context() + call psb_exit(l_ctxt,close=.false.) + end if + call sv%desc_local%free(info) + nullify(sv%a) + call psb_erractionrestore(err_act) + return + end subroutine s_rkr_solver_clear_data + + + subroutine s_rkr_solver_free(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='s_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%clear_data(info) + + call psb_erractionrestore(err_act) + return + end subroutine s_rkr_solver_free + + function s_rkr_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "RKR solver" + end function s_rkr_solver_get_fmt + + subroutine s_rkr_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(amg_s_rkr_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='amg_s_rkr_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = psb_out_unit + endif + + if (sv%global) then + write(iout_,*) ' Recursive Krylov solver (global)' + else + write(iout_,*) ' Recursive Krylov solver (local) ' + end if + write(iout_,*) ' method: ',sv%method + write(iout_,*) ' kprec: ',sv%kprec + if (sv%i_sub_solve > 0) then + write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve) + else + write(iout_,*) ' sub_solve: ',sv%sub_solve + end if + write(iout_,*) ' itmax: ',sv%itmax + write(iout_,*) ' eps: ',sv%eps + write(iout_,*) ' fillin: ',sv%fillin + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine s_rkr_solver_descr + + subroutine s_rkr_solver_cnv(sv,info,amold,vmold,imold) + implicit none + class(amg_s_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + + call sv%prec%cnv(info,amold=amold,vmold=vmold,imold=imold) + + end subroutine s_rkr_solver_cnv + + subroutine s_rkr_solver_clone(sv,svout,info) + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + class(amg_s_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + call svout%free(info) + allocate(svout,stat=info,mold=sv) + select type(so=>svout) + class is(amg_s_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + call sv%a_local%clone(so%a_local,info) + call sv%desc_local%clone(so%desc_local,info) + call sv%prec%clone(so%prec,info) + class default + info = psb_err_internal_error_ + end select + + end subroutine s_rkr_solver_clone + + + subroutine s_rkr_solver_clone_settings(sv,svout,info) + Implicit None + + ! Arguments + class(amg_s_rkr_solver_type), intent(inout) :: sv + class(amg_s_base_solver_type), intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + select type(so=>svout) + class is(amg_s_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + class default + info = psb_err_internal_error_ + end select + + end subroutine s_rkr_solver_clone_settings + + subroutine s_rkr_solver_dmp(sv,desc,level,info,prefix,head,solver,global_num) + implicit none + class(amg_s_rkr_solver_type), intent(in) :: sv + type(psb_desc_type), intent(in) :: desc + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver, global_num + + + call sv%prec%dump(info,prefix=prefix,head=head) + + end subroutine s_rkr_solver_dmp + ! + ! Notify whether RKR is used as a global solver + ! + function s_rkr_solver_is_global(sv) result(val) + implicit none + class(amg_s_rkr_solver_type), intent(in) :: sv + logical :: val + + val = (sv%global) + end function s_rkr_solver_is_global + ! + function s_rkr_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .true. + end function s_rkr_solver_is_iterative + +end module amg_s_rkr_solver diff --git a/amgprec/amg_z_rkr_solver.f90 b/amgprec/amg_z_rkr_solver.f90 new file mode 100644 index 00000000..1eed6c8d --- /dev/null +++ b/amgprec/amg_z_rkr_solver.f90 @@ -0,0 +1,590 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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. +! +! +! +! File: amg_z_rkr_solver_mod.f90 +! +! Module: amg_z_rkr_solver_mod +! +module amg_z_rkr_solver + + use amg_z_base_solver_mod + use amg_z_prec_type + + type, extends(amg_z_base_solver_type) :: amg_z_rkr_solver_type + ! + logical :: global + character(len=16) :: method, kprec, sub_solve + integer(psb_ipk_) :: irst, istopc, itmax, itrace, i_sub_solve + integer(psb_ipk_) :: fillin + real(psb_dpk_) :: eps + type(amg_zprec_type) :: prec + type(psb_desc_type) :: desc_local + type(psb_zspmat_type) :: a_local + type(psb_z_vect_type) :: x_local, z_local + type(psb_zspmat_type), pointer :: a=>null() + contains + ! + ! + procedure, pass(sv) :: dump => z_rkr_solver_dmp + procedure, pass(sv) :: check => z_rkr_solver_check + procedure, pass(sv) :: clone => z_rkr_solver_clone + procedure, pass(sv) :: clone_settings => z_rkr_solver_clone_settings + procedure, pass(sv) :: cnv => z_rkr_solver_cnv + procedure, pass(sv) :: apply_v => amg_z_rkr_solver_apply_vect + procedure, pass(sv) :: apply_a => amg_z_rkr_solver_apply + procedure, pass(sv) :: clear_data => z_rkr_solver_clear_data + procedure, pass(sv) :: free => z_rkr_solver_free + procedure, pass(sv) :: cseti => z_rkr_solver_cseti + procedure, pass(sv) :: csetc => z_rkr_solver_csetc + procedure, pass(sv) :: csetr => z_rkr_solver_csetr + procedure, pass(sv) :: sizeof => z_rkr_solver_sizeof + procedure, pass(sv) :: get_nzeros => z_rkr_solver_get_nzeros + !procedure, nopass :: get_id => z_rkr_solver_get_id + procedure, pass(sv) :: is_global => z_rkr_solver_is_global + procedure, nopass :: is_iterative => z_rkr_solver_is_iterative + + + ! + ! These methods are specific for the new solver type + ! and therefore need to be overridden + ! + procedure, pass(sv) :: descr => z_rkr_solver_descr + procedure, pass(sv) :: default => z_rkr_solver_default + procedure, pass(sv) :: build => amg_z_rkr_solver_bld + procedure, nopass :: get_fmt => z_rkr_solver_get_fmt + end type amg_z_rkr_solver_type + + + private :: z_rkr_solver_get_fmt, z_rkr_solver_descr, z_rkr_solver_default + + interface + subroutine amg_z_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + import :: psb_desc_type, amg_z_rkr_solver_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_rkr_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_z_vect_type),intent(inout), optional :: initu + end subroutine amg_z_rkr_solver_apply_vect + end interface + + interface + subroutine amg_z_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + import :: psb_desc_type, amg_z_rkr_solver_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type, psb_ipk_ + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_rkr_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_dpk_),intent(inout), optional :: initu(:) + end subroutine amg_z_rkr_solver_apply + end interface + + interface + subroutine amg_z_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) + import :: psb_desc_type, amg_z_rkr_solver_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_z_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_z_rkr_solver_bld + end interface + + +contains + + ! + ! + subroutine z_rkr_solver_default(sv) + + Implicit None + + ! Arguments + class(amg_z_rkr_solver_type), intent(inout) :: sv + + sv%method = 'bicgstab' + sv%kprec = 'bjac' + sv%sub_solve = 'ilu' + sv%i_sub_solve = -1 + sv%fillin = 1 + sv%irst = 30 + sv%istopc = 2 + sv%itmax = 40 + sv%itrace = -1 + sv%eps = 1.d-6 + sv%global = .false. + + return + end subroutine z_rkr_solver_default + + function z_rkr_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(amg_z_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%get_nzeros() + + return + end function z_rkr_solver_get_nzeros + + function z_rkr_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(amg_z_rkr_solver_type), intent(in) :: sv + integer(psb_epk_) :: val + + val = sv%prec%sizeof() + sv%desc_local%sizeof() + sv%a_local%sizeof() + + return + end function z_rkr_solver_sizeof + + + subroutine z_rkr_solver_check(sv,info) + + Implicit None + + ! Arguments + class(amg_z_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='z_rkr_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine z_rkr_solver_check + + subroutine z_rkr_solver_cseti(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_z_rkr_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_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='z_rkr_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(trim(what))) + case('RKR_IRST') + sv%irst = val + case('RKR_ISTOPC') + sv%istopc = val + case('RKR_ITMAX') + sv%itmax = val + case('RKR_ITRACE') + sv%itrace = val + case('RKR_SUB_SOLVE') + sv%i_sub_solve = val + case('RKR_FILLIN') + sv%fillin = val + case default + call sv%amg_z_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine z_rkr_solver_cseti + + subroutine z_rkr_solver_csetc(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_z_rkr_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='z_rkr_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + select case(psb_toupper(trim(what))) + case('RKR_METHOD') + sv%method = psb_toupper(trim(val)) + case('RKR_KPREC') + sv%kprec = psb_toupper(trim(val)) + case('RKR_SUB_SOLVE') + sv%sub_solve = psb_toupper(trim(val)) + case('RKR_GLOBAL') + select case(psb_toupper(trim(val))) + case('LOCAL','FALSE') + sv%global = .false. + case('GLOBAL','TRUE') + sv%global =.true. + end select + case default + call sv%amg_z_base_solver_type%set(what,val,info,idx=idx) + end select + + + 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 z_rkr_solver_csetc + + subroutine z_rkr_solver_csetr(sv,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_z_rkr_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_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='z_rkr_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('RKR_EPS') + sv%eps = val + case default + call sv%amg_z_base_solver_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine z_rkr_solver_csetr + + subroutine z_rkr_solver_clear_data(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_z_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='z_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%prec%free(info) + call sv%a_local%free() + call sv%x_local%free(info) + call sv%z_local%free(info) + if ((.not.sv%global).and.sv%desc_local%is_ok()) then + l_ctxt=sv%desc_local%get_context() + call psb_exit(l_ctxt,close=.false.) + end if + call sv%desc_local%free(info) + nullify(sv%a) + call psb_erractionrestore(err_act) + return + end subroutine z_rkr_solver_clear_data + + + subroutine z_rkr_solver_free(sv,info) + use psb_base_mod, only : psb_exit + Implicit None + + ! Arguments + class(amg_z_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: l_ctxt + character(len=20) :: name='z_rkr_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%clear_data(info) + + call psb_erractionrestore(err_act) + return + end subroutine z_rkr_solver_free + + function z_rkr_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "RKR solver" + end function z_rkr_solver_get_fmt + + subroutine z_rkr_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(amg_z_rkr_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='amg_z_rkr_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = psb_out_unit + endif + + if (sv%global) then + write(iout_,*) ' Recursive Krylov solver (global)' + else + write(iout_,*) ' Recursive Krylov solver (local) ' + end if + write(iout_,*) ' method: ',sv%method + write(iout_,*) ' kprec: ',sv%kprec + if (sv%i_sub_solve > 0) then + write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve) + else + write(iout_,*) ' sub_solve: ',sv%sub_solve + end if + write(iout_,*) ' itmax: ',sv%itmax + write(iout_,*) ' eps: ',sv%eps + write(iout_,*) ' fillin: ',sv%fillin + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine z_rkr_solver_descr + + subroutine z_rkr_solver_cnv(sv,info,amold,vmold,imold) + implicit none + class(amg_z_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + + call sv%prec%cnv(info,amold=amold,vmold=vmold,imold=imold) + + end subroutine z_rkr_solver_cnv + + subroutine z_rkr_solver_clone(sv,svout,info) + Implicit None + + ! Arguments + class(amg_z_rkr_solver_type), intent(inout) :: sv + class(amg_z_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + call svout%free(info) + allocate(svout,stat=info,mold=sv) + select type(so=>svout) + class is(amg_z_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + call sv%a_local%clone(so%a_local,info) + call sv%desc_local%clone(so%desc_local,info) + call sv%prec%clone(so%prec,info) + class default + info = psb_err_internal_error_ + end select + + end subroutine z_rkr_solver_clone + + + subroutine z_rkr_solver_clone_settings(sv,svout,info) + Implicit None + + ! Arguments + class(amg_z_rkr_solver_type), intent(inout) :: sv + class(amg_z_base_solver_type), intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + + select type(so=>svout) + class is(amg_z_rkr_solver_type) + so%method = sv%method + so%kprec = sv%kprec + so%sub_solve = sv%sub_solve + so%i_sub_solve = sv%i_sub_solve + so%fillin = sv%fillin + so%irst = sv%irst + so%istopc = sv%istopc + so%itmax = sv%itmax + so%itrace = sv%itrace + so%global = sv%global + so%eps = sv%eps + class default + info = psb_err_internal_error_ + end select + + end subroutine z_rkr_solver_clone_settings + + subroutine z_rkr_solver_dmp(sv,desc,level,info,prefix,head,solver,global_num) + implicit none + class(amg_z_rkr_solver_type), intent(in) :: sv + type(psb_desc_type), intent(in) :: desc + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver, global_num + + + call sv%prec%dump(info,prefix=prefix,head=head) + + end subroutine z_rkr_solver_dmp + ! + ! Notify whether RKR is used as a global solver + ! + function z_rkr_solver_is_global(sv) result(val) + implicit none + class(amg_z_rkr_solver_type), intent(in) :: sv + logical :: val + + val = (sv%global) + end function z_rkr_solver_is_global + ! + function z_rkr_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .true. + end function z_rkr_solver_is_iterative + +end module amg_z_rkr_solver diff --git a/amgprec/impl/solver/Makefile b/amgprec/impl/solver/Makefile index 2ad66c71..497cec4e 100644 --- a/amgprec/impl/solver/Makefile +++ b/amgprec/impl/solver/Makefile @@ -300,7 +300,11 @@ amg_z_invk_solver_check.o \ amg_z_invk_solver_clone.o \ amg_z_invk_solver_cseti.o \ amg_z_invk_solver_descr.o \ -amg_z_invk_solver_seti.o +amg_z_invk_solver_seti.o \ +amg_d_rkr_solver_impl.o \ +amg_s_rkr_solver_impl.o \ +amg_c_rkr_solver_impl.o \ +amg_z_rkr_solver_impl.o LIBNAME=libamg_prec.a diff --git a/amgprec/impl/solver/amg_c_rkr_solver_impl.f90 b/amgprec/impl/solver/amg_c_rkr_solver_impl.f90 new file mode 100644 index 00000000..97ed98cc --- /dev/null +++ b/amgprec/impl/solver/amg_c_rkr_solver_impl.f90 @@ -0,0 +1,271 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! +! File: amg_c_rkr_solver_impl.f90 +! +! This is the implementation file corresponding to amg_c_rkr_solver_mod. +! +! +subroutine amg_c_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold) + + use psb_base_mod + use amg_c_rkr_solver, amg_protect_name => amg_c_rkr_solver_bld + + Implicit None + + ! Arguments + type(psb_cspmat_type), intent(inout), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_c_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_lpk_) :: lnr + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt, l_ctxt + character(len=20) :: name='c_rkr_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call sv%prec%free(info) + + if (sv%global) then + call sv%prec%init(ctxt,sv%kprec,info) + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call sv%prec%hierarchy_build(a,desc_a,info) + call sv%prec%smoothers_build(a,desc_a,info,amold=amold,vmold=vmold) + sv%a => a + else + call psb_init(l_ctxt,np=1_psb_ipk_,basectxt=ctxt,ids=(/me/)) + n_row = desc_a%get_local_rows() + lnr = n_row + call psb_cdall(l_ctxt,sv%desc_local,info,mg=lnr,repl=.true.) + call psb_cdasb(sv%desc_local,info) + call sv%prec%init(l_ctxt,sv%kprec,info) + !This is a gigantic kludge, to be revised + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call a%csclip(sv%a_local,info,jmax=n_row) + call sv%prec%hierarchy_build(sv%a_local,sv%desc_local,info) + call sv%prec%smoothers_build(sv%a_local,sv%desc_local,info,amold=amold,vmold=vmold) + call psb_geall(sv%x_local,sv%desc_local,info) + call sv%x_local%zero() + call psb_geasb(sv%x_local,sv%desc_local,info,mold=vmold) + call psb_geall(sv%z_local,sv%desc_local,info) + call sv%z_local%zero() + call psb_geasb(sv%z_local,sv%desc_local,info,mold=vmold) + end if + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='inner precbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_c_rkr_solver_bld + +subroutine amg_c_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + + use psb_base_mod + use psb_krylov_mod + use amg_c_rkr_solver, amg_protect_name => amg_c_rkr_solver_apply_vect + + Implicit None + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_rkr_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + type(psb_c_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_c_vect_type),intent(inout), optional :: initu + + type(psb_c_vect_type) :: z + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='c_rkr_solver_apply_v', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + if (sv%global) then + call psb_geall(z,desc_data,info) + call z%zero() + call psb_geasb(z,desc_data,info,mold=x%v) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' norms ',& + & psb_genrm2(x,desc_data,info),psb_genrm2(y,desc_data,info),& + & psb_genrm2(z,desc_data,info) + + call psb_krylov(sv%method,sv%a,sv%prec,x,z,sv%eps,& + & desc_data,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) +!!$ call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + else + call psb_geaxpby(cone,x,czero,sv%x_local,sv%desc_local,info) + call psb_krylov(sv%method,sv%a_local,sv%prec,sv%x_local,sv%z_local,sv%eps,& + & sv%desc_local,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) + call psb_geaxpby(alpha,sv%z_local,beta,y,sv%desc_local,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_c_rkr_solver_apply_vect + + +subroutine amg_c_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_c_rkr_solver, amg_protect_name => amg_c_rkr_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_c_rkr_solver_type), intent(inout) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_spk_),intent(inout), optional :: initu(:) + complex(psb_spk_), allocatable :: z(:) + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='c_rkr_solver_apply', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call psb_geasb(z,desc_data,info,scratch=.true.) + + call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_c_rkr_solver_apply diff --git a/amgprec/impl/solver/amg_d_rkr_solver_impl.f90 b/amgprec/impl/solver/amg_d_rkr_solver_impl.f90 new file mode 100644 index 00000000..00adfdde --- /dev/null +++ b/amgprec/impl/solver/amg_d_rkr_solver_impl.f90 @@ -0,0 +1,271 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! +! File: amg_d_rkr_solver_impl.f90 +! +! This is the implementation file corresponding to amg_d_rkr_solver_mod. +! +! +subroutine amg_d_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold) + + use psb_base_mod + use amg_d_rkr_solver, amg_protect_name => amg_d_rkr_solver_bld + + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(inout), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_rkr_solver_type), intent(inout) :: sv + 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 + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_lpk_) :: lnr + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt, l_ctxt + character(len=20) :: name='d_rkr_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call sv%prec%free(info) + + if (sv%global) then + call sv%prec%init(ctxt,sv%kprec,info) + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call sv%prec%hierarchy_build(a,desc_a,info) + call sv%prec%smoothers_build(a,desc_a,info,amold=amold,vmold=vmold) + sv%a => a + else + call psb_init(l_ctxt,np=1_psb_ipk_,basectxt=ctxt,ids=(/me/)) + n_row = desc_a%get_local_rows() + lnr = n_row + call psb_cdall(l_ctxt,sv%desc_local,info,mg=lnr,repl=.true.) + call psb_cdasb(sv%desc_local,info) + call sv%prec%init(l_ctxt,sv%kprec,info) + !This is a gigantic kludge, to be revised + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call a%csclip(sv%a_local,info,jmax=n_row) + call sv%prec%hierarchy_build(sv%a_local,sv%desc_local,info) + call sv%prec%smoothers_build(sv%a_local,sv%desc_local,info,amold=amold,vmold=vmold) + call psb_geall(sv%x_local,sv%desc_local,info) + call sv%x_local%zero() + call psb_geasb(sv%x_local,sv%desc_local,info,mold=vmold) + call psb_geall(sv%z_local,sv%desc_local,info) + call sv%z_local%zero() + call psb_geasb(sv%z_local,sv%desc_local,info,mold=vmold) + end if + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='inner precbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_d_rkr_solver_bld + +subroutine amg_d_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + + use psb_base_mod + use psb_krylov_mod + use amg_d_rkr_solver, amg_protect_name => amg_d_rkr_solver_apply_vect + + Implicit None + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_rkr_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(:) + type(psb_d_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_d_vect_type),intent(inout), optional :: initu + + type(psb_d_vect_type) :: z + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='d_rkr_solver_apply_v', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + if (sv%global) then + call psb_geall(z,desc_data,info) + call z%zero() + call psb_geasb(z,desc_data,info,mold=x%v) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' norms ',& + & psb_genrm2(x,desc_data,info),psb_genrm2(y,desc_data,info),& + & psb_genrm2(z,desc_data,info) + + call psb_krylov(sv%method,sv%a,sv%prec,x,z,sv%eps,& + & desc_data,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) +!!$ call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + else + call psb_geaxpby(done,x,dzero,sv%x_local,sv%desc_local,info) + call psb_krylov(sv%method,sv%a_local,sv%prec,sv%x_local,sv%z_local,sv%eps,& + & sv%desc_local,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) + call psb_geaxpby(alpha,sv%z_local,beta,y,sv%desc_local,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_d_rkr_solver_apply_vect + + +subroutine amg_d_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_d_rkr_solver, amg_protect_name => amg_d_rkr_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_rkr_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 + character, intent(in), optional :: init + real(psb_dpk_),intent(inout), optional :: initu(:) + real(psb_dpk_), allocatable :: z(:) + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='d_rkr_solver_apply', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call psb_geasb(z,desc_data,info,scratch=.true.) + + call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_d_rkr_solver_apply diff --git a/amgprec/impl/solver/amg_s_rkr_solver_impl.f90 b/amgprec/impl/solver/amg_s_rkr_solver_impl.f90 new file mode 100644 index 00000000..492fc088 --- /dev/null +++ b/amgprec/impl/solver/amg_s_rkr_solver_impl.f90 @@ -0,0 +1,271 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! +! File: amg_s_rkr_solver_impl.f90 +! +! This is the implementation file corresponding to amg_s_rkr_solver_mod. +! +! +subroutine amg_s_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold) + + use psb_base_mod + use amg_s_rkr_solver, amg_protect_name => amg_s_rkr_solver_bld + + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(inout), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_s_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_lpk_) :: lnr + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt, l_ctxt + character(len=20) :: name='s_rkr_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call sv%prec%free(info) + + if (sv%global) then + call sv%prec%init(ctxt,sv%kprec,info) + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call sv%prec%hierarchy_build(a,desc_a,info) + call sv%prec%smoothers_build(a,desc_a,info,amold=amold,vmold=vmold) + sv%a => a + else + call psb_init(l_ctxt,np=1_psb_ipk_,basectxt=ctxt,ids=(/me/)) + n_row = desc_a%get_local_rows() + lnr = n_row + call psb_cdall(l_ctxt,sv%desc_local,info,mg=lnr,repl=.true.) + call psb_cdasb(sv%desc_local,info) + call sv%prec%init(l_ctxt,sv%kprec,info) + !This is a gigantic kludge, to be revised + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call a%csclip(sv%a_local,info,jmax=n_row) + call sv%prec%hierarchy_build(sv%a_local,sv%desc_local,info) + call sv%prec%smoothers_build(sv%a_local,sv%desc_local,info,amold=amold,vmold=vmold) + call psb_geall(sv%x_local,sv%desc_local,info) + call sv%x_local%zero() + call psb_geasb(sv%x_local,sv%desc_local,info,mold=vmold) + call psb_geall(sv%z_local,sv%desc_local,info) + call sv%z_local%zero() + call psb_geasb(sv%z_local,sv%desc_local,info,mold=vmold) + end if + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='inner precbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_rkr_solver_bld + +subroutine amg_s_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + + use psb_base_mod + use psb_krylov_mod + use amg_s_rkr_solver, amg_protect_name => amg_s_rkr_solver_apply_vect + + Implicit None + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_rkr_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_s_vect_type),intent(inout), optional :: initu + + type(psb_s_vect_type) :: z + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='s_rkr_solver_apply_v', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + if (sv%global) then + call psb_geall(z,desc_data,info) + call z%zero() + call psb_geasb(z,desc_data,info,mold=x%v) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' norms ',& + & psb_genrm2(x,desc_data,info),psb_genrm2(y,desc_data,info),& + & psb_genrm2(z,desc_data,info) + + call psb_krylov(sv%method,sv%a,sv%prec,x,z,sv%eps,& + & desc_data,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) +!!$ call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + else + call psb_geaxpby(sone,x,szero,sv%x_local,sv%desc_local,info) + call psb_krylov(sv%method,sv%a_local,sv%prec,sv%x_local,sv%z_local,sv%eps,& + & sv%desc_local,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) + call psb_geaxpby(alpha,sv%z_local,beta,y,sv%desc_local,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_rkr_solver_apply_vect + + +subroutine amg_s_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_s_rkr_solver, amg_protect_name => amg_s_rkr_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_rkr_solver_type), intent(inout) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + real(psb_spk_),intent(inout), optional :: initu(:) + real(psb_spk_), allocatable :: z(:) + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='s_rkr_solver_apply', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call psb_geasb(z,desc_data,info,scratch=.true.) + + call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_rkr_solver_apply diff --git a/amgprec/impl/solver/amg_z_rkr_solver_impl.f90 b/amgprec/impl/solver/amg_z_rkr_solver_impl.f90 new file mode 100644 index 00000000..c6d96249 --- /dev/null +++ b/amgprec/impl/solver/amg_z_rkr_solver_impl.f90 @@ -0,0 +1,271 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2020 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! moved here from amg4psblas-extension +! +! +! AMG4PSBLAS Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, IT +! +! 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 AMG4PSBLAS 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 AMG4PSBLAS 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. +! +! +! File: amg_z_rkr_solver_impl.f90 +! +! This is the implementation file corresponding to amg_z_rkr_solver_mod. +! +! +subroutine amg_z_rkr_solver_bld(a,desc_a,sv,info,b,amold,vmold) + + use psb_base_mod + use amg_z_rkr_solver, amg_protect_name => amg_z_rkr_solver_bld + + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(inout), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_z_rkr_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_lpk_) :: lnr + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt, l_ctxt + character(len=20) :: name='z_rkr_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call sv%prec%free(info) + + if (sv%global) then + call sv%prec%init(ctxt,sv%kprec,info) + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call sv%prec%hierarchy_build(a,desc_a,info) + call sv%prec%smoothers_build(a,desc_a,info,amold=amold,vmold=vmold) + sv%a => a + else + call psb_init(l_ctxt,np=1_psb_ipk_,basectxt=ctxt,ids=(/me/)) + n_row = desc_a%get_local_rows() + lnr = n_row + call psb_cdall(l_ctxt,sv%desc_local,info,mg=lnr,repl=.true.) + call psb_cdasb(sv%desc_local,info) + call sv%prec%init(l_ctxt,sv%kprec,info) + !This is a gigantic kludge, to be revised + if (sv%i_sub_solve>0) then + call sv%prec%set('sub_solve',sv%i_sub_solve,info) + else + call sv%prec%set('sub_solve',sv%sub_solve,info) + end if + call sv%prec%set('sub_fillin',sv%fillin,info) + call a%csclip(sv%a_local,info,jmax=n_row) + call sv%prec%hierarchy_build(sv%a_local,sv%desc_local,info) + call sv%prec%smoothers_build(sv%a_local,sv%desc_local,info,amold=amold,vmold=vmold) + call psb_geall(sv%x_local,sv%desc_local,info) + call sv%x_local%zero() + call psb_geasb(sv%x_local,sv%desc_local,info,mold=vmold) + call psb_geall(sv%z_local,sv%desc_local,info) + call sv%z_local%zero() + call psb_geasb(sv%z_local,sv%desc_local,info,mold=vmold) + end if + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='inner precbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_z_rkr_solver_bld + +subroutine amg_z_rkr_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& + & trans,work,wv,info,init,initu) + + use psb_base_mod + use psb_krylov_mod + use amg_z_rkr_solver, amg_protect_name => amg_z_rkr_solver_apply_vect + + Implicit None + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_rkr_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + type(psb_z_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_z_vect_type),intent(inout), optional :: initu + + type(psb_z_vect_type) :: z + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='z_rkr_solver_apply_v', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + if (sv%global) then + call psb_geall(z,desc_data,info) + call z%zero() + call psb_geasb(z,desc_data,info,mold=x%v) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' norms ',& + & psb_genrm2(x,desc_data,info),psb_genrm2(y,desc_data,info),& + & psb_genrm2(z,desc_data,info) + + call psb_krylov(sv%method,sv%a,sv%prec,x,z,sv%eps,& + & desc_data,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) +!!$ call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + else + call psb_geaxpby(zone,x,zzero,sv%x_local,sv%desc_local,info) + call psb_krylov(sv%method,sv%a_local,sv%prec,sv%x_local,sv%z_local,sv%eps,& + & sv%desc_local,info,itmax=sv%itmax,itrace=sv%itrace,& + & istop=sv%istopc,irst=sv%irst) + call psb_geaxpby(alpha,sv%z_local,beta,y,sv%desc_local,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_z_rkr_solver_apply_vect + + +subroutine amg_z_rkr_solver_apply(alpha,sv,x,beta,y,desc_data,& + & trans,work,info,init,initu) + use psb_base_mod + use amg_z_rkr_solver, amg_protect_name => amg_z_rkr_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_z_rkr_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + complex(psb_dpk_),intent(inout), optional :: initu(:) + complex(psb_dpk_), allocatable :: z(:) + integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level + type(psb_ctxt_type) :: ctxt + character(len=20) :: name='z_rkr_solver_apply', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_data%get_context() + call psb_info(ctxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + call psb_geasb(z,desc_data,info,scratch=.true.) + + call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work) + + call psb_geaxpby(alpha,z,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_z_rkr_solver_apply diff --git a/tests/pdegen/amg_d_pde2d.f90 b/tests/pdegen/amg_d_pde2d.f90 index 71630e5b..72c76452 100644 --- a/tests/pdegen/amg_d_pde2d.f90 +++ b/tests/pdegen/amg_d_pde2d.f90 @@ -75,6 +75,8 @@ program amg_d_pde2d use amg_d_genpde_mod use amg_ainv_mod use amg_d_ilu_solver + use amg_d_rkr_solver + implicit none ! input parameters @@ -117,6 +119,7 @@ program amg_d_pde2d type(amg_d_invt_solver_type) :: invtsv type(amg_d_invk_solver_type) :: invksv type(amg_d_ainv_solver_type) :: ainvsv + type(amg_d_rkr_solver_type) :: rkr_slv type precdata ! preconditioner type @@ -173,10 +176,22 @@ program amg_d_pde2d ! (repl. mat.) character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, ! MILU, UMF, MUMPS, SLU + character(len=16) :: cvariant ! AINV variant: LLK, etc + character(len=16) :: ckryl ! Krylov method for RKR, ignored otherwise + ! MILU, UMF, MUMPS, SLU integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization - real(psb_dpk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: cinvfill ! Inverse fill-in for INVK integer(psb_ipk_) :: cjswp ! sweeps for GS or JAC coarsest-lev subsolver - + real(psb_dpk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: crkiter ! Max iterations for RKR, ignored otherwise + real(psb_dpk_) :: crkeps ! eps for RKR, ignored otherwise + integer(psb_ipk_) :: crktrace ! ITRACE for RKR, ignored otherwise + character(len=16) :: checkres ! Check the BJAC residual + character(len=16) :: printres ! Print the BJAC residual + integer(psb_ipk_) :: checkiter ! ITRACE for residual check + integer(psb_ipk_) :: printiter ! ITRACE for residual print + real(psb_dpk_) :: tol ! Tolerance for exit from BJAC + logical :: dump ! dump on file end type precdata type(precdata) :: p_choice @@ -354,13 +369,72 @@ program amg_d_pde2d end select end if - call prec%set('coarse_solve', p_choice%csolve, info) - if (psb_toupper(p_choice%csolve) == 'BJAC') & - & call prec%set('coarse_subsolve', p_choice%csbsolve, info) - call prec%set('coarse_mat', p_choice%cmat, info) - call prec%set('coarse_fillin', p_choice%cfill, info) - call prec%set('coarse_iluthrs', p_choice%cthres, info) - call prec%set('coarse_sweeps', p_choice%cjswp, info) + nlv = size(prec%precv) + if (psb_toupper(p_choice%csolve) == 'RKR') then + call prec%set('coarse_solve', 'BJAC', info) + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_kprec', 'BJAC' , info,ilev=nlv) + call prec%set('rkr_global', 'global' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + else + call prec%set('coarse_solve', p_choice%csolve, info) + select case(psb_toupper(trim(p_choice%csolve))) + case('BJAC','L1-BJAC') + select case(trim(psb_toupper(p_choice%csbsolve))) + case ('MUMPS') + call prec%set('MUMPS_LOC_GLOB','LOCAL_SOLVER',info) + ! + ! Experimental: enable Block Low-Rank + ! + call prec%set('MUMPS_IPAR_ENTRY',ione,info,idx=35_psb_ipk_) + call prec%set('MUMPS_RPAR_ENTRY',4.0_psb_dpk_,info,idx=7_psb_ipk_) + case('RKR') + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_global', 'local' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + case('INVK') + call prec%set(invksv, info,ilev=nlv) + case('INVT') + call prec%set(invtsv, info,ilev=nlv) + case('AINV') + call prec%set(ainvsv, info,ilev=nlv) + call prec%set('ainv_alg', p_choice%cvariant, info, ilev=nlv) + case default + call prec%set('coarse_subsolve', p_choice%csbsolve, info) + end select + ! + ! Experimental: Use residual check on coarse BJAC solver + ! + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case('L1-GS') + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case default + ! Do nothing, should have been already handled + end select + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_fillin', p_choice%cfill, info) + call prec%set('inv_fillin', p_choice%cinvfill, info, ilev=nlv) + call prec%set('coarse_iluthrs', p_choice%cthres, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + end if end select @@ -574,10 +648,21 @@ contains ! coasest-level solver call read_data(prec%csolve,inp_unit) ! coarsest-lev solver call read_data(prec%csbsolve,inp_unit) ! coarsest-lev subsolver + call read_data(prec%cvariant,inp_unit) ! AINV variant + call read_data(prec%ckryl,inp_unit) ! RKR Krylov method choice call read_data(prec%cmat,inp_unit) ! coarsest mat layout call read_data(prec%cfill,inp_unit) ! fill-in for incompl LU + call read_data(prec%cinvfill,inp_unit) ! Inverse fill-in for INVK call read_data(prec%cthres,inp_unit) ! Threshold for ILUT call read_data(prec%cjswp,inp_unit) ! sweeps for GS/JAC subsolver + call read_data(prec%crkiter,inp_unit) ! Max iterations for RKR, ignored otherwise + call read_data(prec%crkeps,inp_unit) ! eps for RKR, ignored otherwise + call read_data(prec%crktrace,inp_unit) ! itrace for RKR, ignored otherwise + call read_data(prec%checkres,inp_unit) ! Check the BJAC residual + call read_data(prec%printres,inp_unit) ! Print the BJAC residual + call read_data(prec%checkiter,inp_unit) ! ITRACE for residual check + call read_data(prec%printiter,inp_unit) ! ITRACE for residual print + call read_data(prec%tol,inp_unit) ! Tolerance for exit from BJAC if (inp_unit /= psb_inp_unit) then close(inp_unit) end if @@ -637,14 +722,26 @@ contains end if call psb_bcast(ctxt,prec%athres) + ! broadcast coasest-level solver call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) call psb_bcast(ctxt,prec%csbsolve) + call psb_bcast(ctxt,prec%cvariant) + call psb_bcast(ctxt,prec%ckryl) call psb_bcast(ctxt,prec%cfill) + call psb_bcast(ctxt,prec%cinvfill) call psb_bcast(ctxt,prec%cthres) call psb_bcast(ctxt,prec%cjswp) - + call psb_bcast(ctxt,prec%crkiter) + call psb_bcast(ctxt,prec%crkeps) + call psb_bcast(ctxt,prec%crktrace) + call psb_bcast(ctxt,prec%checkres) + call psb_bcast(ctxt,prec%printres) + call psb_bcast(ctxt,prec%checkiter) + call psb_bcast(ctxt,prec%printiter) + call psb_bcast(ctxt,prec%tol) + call psb_bcast(ctxt,prec%dump) end subroutine get_parms diff --git a/tests/pdegen/amg_d_pde3d.f90 b/tests/pdegen/amg_d_pde3d.f90 index 13323797..a2d2139e 100644 --- a/tests/pdegen/amg_d_pde3d.f90 +++ b/tests/pdegen/amg_d_pde3d.f90 @@ -76,6 +76,8 @@ program amg_d_pde3d use amg_d_genpde_mod use amg_ainv_mod use amg_d_ilu_solver + use amg_d_rkr_solver + implicit none ! input parameters @@ -118,6 +120,7 @@ program amg_d_pde3d type(amg_d_invt_solver_type) :: invtsv type(amg_d_invk_solver_type) :: invksv type(amg_d_ainv_solver_type) :: ainvsv + type(amg_d_rkr_solver_type) :: rkr_slv type precdata ! preconditioner type @@ -174,10 +177,22 @@ program amg_d_pde3d ! (repl. mat.) character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, ! MILU, UMF, MUMPS, SLU + character(len=16) :: cvariant ! AINV variant: LLK, etc + character(len=16) :: ckryl ! Krylov method for RKR, ignored otherwise + ! MILU, UMF, MUMPS, SLU integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization - real(psb_dpk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: cinvfill ! Inverse fill-in for INVK integer(psb_ipk_) :: cjswp ! sweeps for GS or JAC coarsest-lev subsolver - + real(psb_dpk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: crkiter ! Max iterations for RKR, ignored otherwise + real(psb_dpk_) :: crkeps ! eps for RKR, ignored otherwise + integer(psb_ipk_) :: crktrace ! ITRACE for RKR, ignored otherwise + character(len=16) :: checkres ! Check the BJAC residual + character(len=16) :: printres ! Print the BJAC residual + integer(psb_ipk_) :: checkiter ! ITRACE for residual check + integer(psb_ipk_) :: printiter ! ITRACE for residual print + real(psb_dpk_) :: tol ! Tolerance for exit from BJAC + logical :: dump ! dump on file end type precdata type(precdata) :: p_choice @@ -358,13 +373,72 @@ program amg_d_pde3d end select end if - call prec%set('coarse_solve', p_choice%csolve, info) - if (psb_toupper(p_choice%csolve) == 'BJAC') & - & call prec%set('coarse_subsolve', p_choice%csbsolve, info) - call prec%set('coarse_mat', p_choice%cmat, info) - call prec%set('coarse_fillin', p_choice%cfill, info) - call prec%set('coarse_iluthrs', p_choice%cthres, info) - call prec%set('coarse_sweeps', p_choice%cjswp, info) + nlv = size(prec%precv) + if (psb_toupper(p_choice%csolve) == 'RKR') then + call prec%set('coarse_solve', 'BJAC', info) + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_kprec', 'BJAC' , info,ilev=nlv) + call prec%set('rkr_global', 'global' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + else + call prec%set('coarse_solve', p_choice%csolve, info) + select case(psb_toupper(trim(p_choice%csolve))) + case('BJAC','L1-BJAC') + select case(trim(psb_toupper(p_choice%csbsolve))) + case ('MUMPS') + call prec%set('MUMPS_LOC_GLOB','LOCAL_SOLVER',info) + ! + ! Experimental: enable Block Low-Rank + ! + call prec%set('MUMPS_IPAR_ENTRY',ione,info,idx=35_psb_ipk_) + call prec%set('MUMPS_RPAR_ENTRY',4.0_psb_dpk_,info,idx=7_psb_ipk_) + case('RKR') + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_global', 'local' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + case('INVK') + call prec%set(invksv, info,ilev=nlv) + case('INVT') + call prec%set(invtsv, info,ilev=nlv) + case('AINV') + call prec%set(ainvsv, info,ilev=nlv) + call prec%set('ainv_alg', p_choice%cvariant, info, ilev=nlv) + case default + call prec%set('coarse_subsolve', p_choice%csbsolve, info) + end select + ! + ! Experimental: Use residual check on coarse BJAC solver + ! + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case('L1-GS') + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case default + ! Do nothing, should have been already handled + end select + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_fillin', p_choice%cfill, info) + call prec%set('inv_fillin', p_choice%cinvfill, info, ilev=nlv) + call prec%set('coarse_iluthrs', p_choice%cthres, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + end if end select @@ -578,10 +652,23 @@ contains ! coasest-level solver call read_data(prec%csolve,inp_unit) ! coarsest-lev solver call read_data(prec%csbsolve,inp_unit) ! coarsest-lev subsolver + call read_data(prec%cvariant,inp_unit) ! AINV variant + call read_data(prec%ckryl,inp_unit) ! RKR Krylov method choice call read_data(prec%cmat,inp_unit) ! coarsest mat layout call read_data(prec%cfill,inp_unit) ! fill-in for incompl LU + call read_data(prec%cinvfill,inp_unit) ! Inverse fill-in for INVK call read_data(prec%cthres,inp_unit) ! Threshold for ILUT call read_data(prec%cjswp,inp_unit) ! sweeps for GS/JAC subsolver + call read_data(prec%crkiter,inp_unit) ! Max iterations for RKR, ignored otherwise + call read_data(prec%crkeps,inp_unit) ! eps for RKR, ignored otherwise + call read_data(prec%crktrace,inp_unit) ! itrace for RKR, ignored otherwise + call read_data(prec%checkres,inp_unit) ! Check the BJAC residual + call read_data(prec%printres,inp_unit) ! Print the BJAC residual + call read_data(prec%checkiter,inp_unit) ! ITRACE for residual check + call read_data(prec%printiter,inp_unit) ! ITRACE for residual print + call read_data(prec%tol,inp_unit) ! Tolerance for exit from BJAC + + call read_data(prec%dump,inp_unit) ! if (inp_unit /= psb_inp_unit) then close(inp_unit) end if @@ -641,13 +728,26 @@ contains end if call psb_bcast(ctxt,prec%athres) + ! broadcast coasest-level solver call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) call psb_bcast(ctxt,prec%csbsolve) + call psb_bcast(ctxt,prec%cvariant) + call psb_bcast(ctxt,prec%ckryl) call psb_bcast(ctxt,prec%cfill) + call psb_bcast(ctxt,prec%cinvfill) call psb_bcast(ctxt,prec%cthres) call psb_bcast(ctxt,prec%cjswp) + call psb_bcast(ctxt,prec%crkiter) + call psb_bcast(ctxt,prec%crkeps) + call psb_bcast(ctxt,prec%crktrace) + call psb_bcast(ctxt,prec%checkres) + call psb_bcast(ctxt,prec%printres) + call psb_bcast(ctxt,prec%checkiter) + call psb_bcast(ctxt,prec%printiter) + call psb_bcast(ctxt,prec%tol) + call psb_bcast(ctxt,prec%dump) end subroutine get_parms diff --git a/tests/pdegen/amg_s_pde2d.f90 b/tests/pdegen/amg_s_pde2d.f90 index 4d380e26..48543a79 100644 --- a/tests/pdegen/amg_s_pde2d.f90 +++ b/tests/pdegen/amg_s_pde2d.f90 @@ -75,6 +75,8 @@ program amg_s_pde2d use amg_s_genpde_mod use amg_ainv_mod use amg_s_ilu_solver + use amg_s_rkr_solver + implicit none ! input parameters @@ -117,6 +119,7 @@ program amg_s_pde2d type(amg_s_invt_solver_type) :: invtsv type(amg_s_invk_solver_type) :: invksv type(amg_s_ainv_solver_type) :: ainvsv + type(amg_s_rkr_solver_type) :: rkr_slv type precdata ! preconditioner type @@ -173,10 +176,22 @@ program amg_s_pde2d ! (repl. mat.) character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, ! MILU, UMF, MUMPS, SLU + character(len=16) :: cvariant ! AINV variant: LLK, etc + character(len=16) :: ckryl ! Krylov method for RKR, ignored otherwise + ! MILU, UMF, MUMPS, SLU integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization - real(psb_spk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: cinvfill ! Inverse fill-in for INVK integer(psb_ipk_) :: cjswp ! sweeps for GS or JAC coarsest-lev subsolver - + real(psb_spk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: crkiter ! Max iterations for RKR, ignored otherwise + real(psb_spk_) :: crkeps ! eps for RKR, ignored otherwise + integer(psb_ipk_) :: crktrace ! ITRACE for RKR, ignored otherwise + character(len=16) :: checkres ! Check the BJAC residual + character(len=16) :: printres ! Print the BJAC residual + integer(psb_ipk_) :: checkiter ! ITRACE for residual check + integer(psb_ipk_) :: printiter ! ITRACE for residual print + real(psb_spk_) :: tol ! Tolerance for exit from BJAC + logical :: dump ! dump on file end type precdata type(precdata) :: p_choice @@ -354,13 +369,72 @@ program amg_s_pde2d end select end if - call prec%set('coarse_solve', p_choice%csolve, info) - if (psb_toupper(p_choice%csolve) == 'BJAC') & - & call prec%set('coarse_subsolve', p_choice%csbsolve, info) - call prec%set('coarse_mat', p_choice%cmat, info) - call prec%set('coarse_fillin', p_choice%cfill, info) - call prec%set('coarse_iluthrs', p_choice%cthres, info) - call prec%set('coarse_sweeps', p_choice%cjswp, info) + nlv = size(prec%precv) + if (psb_toupper(p_choice%csolve) == 'RKR') then + call prec%set('coarse_solve', 'BJAC', info) + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_kprec', 'BJAC' , info,ilev=nlv) + call prec%set('rkr_global', 'global' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + else + call prec%set('coarse_solve', p_choice%csolve, info) + select case(psb_toupper(trim(p_choice%csolve))) + case('BJAC','L1-BJAC') + select case(trim(psb_toupper(p_choice%csbsolve))) + case ('MUMPS') + call prec%set('MUMPS_LOC_GLOB','LOCAL_SOLVER',info) + ! + ! Experimental: enable Block Low-Rank + ! + call prec%set('MUMPS_IPAR_ENTRY',ione,info,idx=35_psb_ipk_) + call prec%set('MUMPS_RPAR_ENTRY',4.0_psb_spk_,info,idx=7_psb_ipk_) + case('RKR') + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_global', 'local' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + case('INVK') + call prec%set(invksv, info,ilev=nlv) + case('INVT') + call prec%set(invtsv, info,ilev=nlv) + case('AINV') + call prec%set(ainvsv, info,ilev=nlv) + call prec%set('ainv_alg', p_choice%cvariant, info, ilev=nlv) + case default + call prec%set('coarse_subsolve', p_choice%csbsolve, info) + end select + ! + ! Experimental: Use residual check on coarse BJAC solver + ! + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case('L1-GS') + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case default + ! Do nothing, should have been already handled + end select + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_fillin', p_choice%cfill, info) + call prec%set('inv_fillin', p_choice%cinvfill, info, ilev=nlv) + call prec%set('coarse_iluthrs', p_choice%cthres, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + end if end select @@ -574,10 +648,21 @@ contains ! coasest-level solver call read_data(prec%csolve,inp_unit) ! coarsest-lev solver call read_data(prec%csbsolve,inp_unit) ! coarsest-lev subsolver + call read_data(prec%cvariant,inp_unit) ! AINV variant + call read_data(prec%ckryl,inp_unit) ! RKR Krylov method choice call read_data(prec%cmat,inp_unit) ! coarsest mat layout call read_data(prec%cfill,inp_unit) ! fill-in for incompl LU + call read_data(prec%cinvfill,inp_unit) ! Inverse fill-in for INVK call read_data(prec%cthres,inp_unit) ! Threshold for ILUT call read_data(prec%cjswp,inp_unit) ! sweeps for GS/JAC subsolver + call read_data(prec%crkiter,inp_unit) ! Max iterations for RKR, ignored otherwise + call read_data(prec%crkeps,inp_unit) ! eps for RKR, ignored otherwise + call read_data(prec%crktrace,inp_unit) ! itrace for RKR, ignored otherwise + call read_data(prec%checkres,inp_unit) ! Check the BJAC residual + call read_data(prec%printres,inp_unit) ! Print the BJAC residual + call read_data(prec%checkiter,inp_unit) ! ITRACE for residual check + call read_data(prec%printiter,inp_unit) ! ITRACE for residual print + call read_data(prec%tol,inp_unit) ! Tolerance for exit from BJAC if (inp_unit /= psb_inp_unit) then close(inp_unit) end if @@ -637,14 +722,26 @@ contains end if call psb_bcast(ctxt,prec%athres) + ! broadcast coasest-level solver call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) call psb_bcast(ctxt,prec%csbsolve) + call psb_bcast(ctxt,prec%cvariant) + call psb_bcast(ctxt,prec%ckryl) call psb_bcast(ctxt,prec%cfill) + call psb_bcast(ctxt,prec%cinvfill) call psb_bcast(ctxt,prec%cthres) call psb_bcast(ctxt,prec%cjswp) - + call psb_bcast(ctxt,prec%crkiter) + call psb_bcast(ctxt,prec%crkeps) + call psb_bcast(ctxt,prec%crktrace) + call psb_bcast(ctxt,prec%checkres) + call psb_bcast(ctxt,prec%printres) + call psb_bcast(ctxt,prec%checkiter) + call psb_bcast(ctxt,prec%printiter) + call psb_bcast(ctxt,prec%tol) + call psb_bcast(ctxt,prec%dump) end subroutine get_parms diff --git a/tests/pdegen/amg_s_pde3d.f90 b/tests/pdegen/amg_s_pde3d.f90 index c17c29db..186f7fc7 100644 --- a/tests/pdegen/amg_s_pde3d.f90 +++ b/tests/pdegen/amg_s_pde3d.f90 @@ -76,6 +76,8 @@ program amg_s_pde3d use amg_s_genpde_mod use amg_ainv_mod use amg_s_ilu_solver + use amg_s_rkr_solver + implicit none ! input parameters @@ -118,6 +120,7 @@ program amg_s_pde3d type(amg_s_invt_solver_type) :: invtsv type(amg_s_invk_solver_type) :: invksv type(amg_s_ainv_solver_type) :: ainvsv + type(amg_s_rkr_solver_type) :: rkr_slv type precdata ! preconditioner type @@ -174,10 +177,22 @@ program amg_s_pde3d ! (repl. mat.) character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, ! MILU, UMF, MUMPS, SLU + character(len=16) :: cvariant ! AINV variant: LLK, etc + character(len=16) :: ckryl ! Krylov method for RKR, ignored otherwise + ! MILU, UMF, MUMPS, SLU integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization - real(psb_spk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: cinvfill ! Inverse fill-in for INVK integer(psb_ipk_) :: cjswp ! sweeps for GS or JAC coarsest-lev subsolver - + real(psb_spk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: crkiter ! Max iterations for RKR, ignored otherwise + real(psb_spk_) :: crkeps ! eps for RKR, ignored otherwise + integer(psb_ipk_) :: crktrace ! ITRACE for RKR, ignored otherwise + character(len=16) :: checkres ! Check the BJAC residual + character(len=16) :: printres ! Print the BJAC residual + integer(psb_ipk_) :: checkiter ! ITRACE for residual check + integer(psb_ipk_) :: printiter ! ITRACE for residual print + real(psb_spk_) :: tol ! Tolerance for exit from BJAC + logical :: dump ! dump on file end type precdata type(precdata) :: p_choice @@ -358,13 +373,72 @@ program amg_s_pde3d end select end if - call prec%set('coarse_solve', p_choice%csolve, info) - if (psb_toupper(p_choice%csolve) == 'BJAC') & - & call prec%set('coarse_subsolve', p_choice%csbsolve, info) - call prec%set('coarse_mat', p_choice%cmat, info) - call prec%set('coarse_fillin', p_choice%cfill, info) - call prec%set('coarse_iluthrs', p_choice%cthres, info) - call prec%set('coarse_sweeps', p_choice%cjswp, info) + nlv = size(prec%precv) + if (psb_toupper(p_choice%csolve) == 'RKR') then + call prec%set('coarse_solve', 'BJAC', info) + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_kprec', 'BJAC' , info,ilev=nlv) + call prec%set('rkr_global', 'global' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + else + call prec%set('coarse_solve', p_choice%csolve, info) + select case(psb_toupper(trim(p_choice%csolve))) + case('BJAC','L1-BJAC') + select case(trim(psb_toupper(p_choice%csbsolve))) + case ('MUMPS') + call prec%set('MUMPS_LOC_GLOB','LOCAL_SOLVER',info) + ! + ! Experimental: enable Block Low-Rank + ! + call prec%set('MUMPS_IPAR_ENTRY',ione,info,idx=35_psb_ipk_) + call prec%set('MUMPS_RPAR_ENTRY',4.0_psb_spk_,info,idx=7_psb_ipk_) + case('RKR') + call prec%set(rkr_slv,info,ilev=nlv) + call prec%set('rkr_method',p_choice%ckryl,info,ilev=nlv) + call prec%set('rkr_global', 'local' , info,ilev=nlv) + call prec%set('rkr_eps', p_choice%crkeps,info,ilev=nlv) + call prec%set('rkr_itmax', p_choice%crkiter,info,ilev=nlv) + call prec%set('rkr_itrace', p_choice%crktrace, info,ilev=nlv) + call prec%set('rkr_fillin', p_choice%cfill,info,ilev=nlv) + case('INVK') + call prec%set(invksv, info,ilev=nlv) + case('INVT') + call prec%set(invtsv, info,ilev=nlv) + case('AINV') + call prec%set(ainvsv, info,ilev=nlv) + call prec%set('ainv_alg', p_choice%cvariant, info, ilev=nlv) + case default + call prec%set('coarse_subsolve', p_choice%csbsolve, info) + end select + ! + ! Experimental: Use residual check on coarse BJAC solver + ! + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case('L1-GS') + call prec%set('SMOOTHER_STOP', p_choice%checkres, info) + call prec%set('SMOOTHER_TRACE', p_choice%printres, info) + call prec%set('SMOOTHER_STOPTOL', p_choice%tol, info) + call prec%set('SMOOTHER_ITRACE', p_choice%printiter, info) + call prec%set('SMOOTHER_RESIDUAL',p_choice%checkiter, info) + case default + ! Do nothing, should have been already handled + end select + call prec%set('coarse_mat', p_choice%cmat, info) + call prec%set('coarse_fillin', p_choice%cfill, info) + call prec%set('inv_fillin', p_choice%cinvfill, info, ilev=nlv) + call prec%set('coarse_iluthrs', p_choice%cthres, info) + call prec%set('coarse_sweeps', p_choice%cjswp, info) + end if end select @@ -578,10 +652,23 @@ contains ! coasest-level solver call read_data(prec%csolve,inp_unit) ! coarsest-lev solver call read_data(prec%csbsolve,inp_unit) ! coarsest-lev subsolver + call read_data(prec%cvariant,inp_unit) ! AINV variant + call read_data(prec%ckryl,inp_unit) ! RKR Krylov method choice call read_data(prec%cmat,inp_unit) ! coarsest mat layout call read_data(prec%cfill,inp_unit) ! fill-in for incompl LU + call read_data(prec%cinvfill,inp_unit) ! Inverse fill-in for INVK call read_data(prec%cthres,inp_unit) ! Threshold for ILUT call read_data(prec%cjswp,inp_unit) ! sweeps for GS/JAC subsolver + call read_data(prec%crkiter,inp_unit) ! Max iterations for RKR, ignored otherwise + call read_data(prec%crkeps,inp_unit) ! eps for RKR, ignored otherwise + call read_data(prec%crktrace,inp_unit) ! itrace for RKR, ignored otherwise + call read_data(prec%checkres,inp_unit) ! Check the BJAC residual + call read_data(prec%printres,inp_unit) ! Print the BJAC residual + call read_data(prec%checkiter,inp_unit) ! ITRACE for residual check + call read_data(prec%printiter,inp_unit) ! ITRACE for residual print + call read_data(prec%tol,inp_unit) ! Tolerance for exit from BJAC + + call read_data(prec%dump,inp_unit) ! if (inp_unit /= psb_inp_unit) then close(inp_unit) end if @@ -641,13 +728,26 @@ contains end if call psb_bcast(ctxt,prec%athres) + ! broadcast coasest-level solver call psb_bcast(ctxt,prec%csize) call psb_bcast(ctxt,prec%cmat) call psb_bcast(ctxt,prec%csolve) call psb_bcast(ctxt,prec%csbsolve) + call psb_bcast(ctxt,prec%cvariant) + call psb_bcast(ctxt,prec%ckryl) call psb_bcast(ctxt,prec%cfill) + call psb_bcast(ctxt,prec%cinvfill) call psb_bcast(ctxt,prec%cthres) call psb_bcast(ctxt,prec%cjswp) + call psb_bcast(ctxt,prec%crkiter) + call psb_bcast(ctxt,prec%crkeps) + call psb_bcast(ctxt,prec%crktrace) + call psb_bcast(ctxt,prec%checkres) + call psb_bcast(ctxt,prec%printres) + call psb_bcast(ctxt,prec%checkiter) + call psb_bcast(ctxt,prec%printiter) + call psb_bcast(ctxt,prec%tol) + call psb_bcast(ctxt,prec%dump) end subroutine get_parms diff --git a/tests/pdegen/runs/amg_pde2d.inp b/tests/pdegen/runs/amg_pde2d.inp index 61d0bcec..07a9901d 100644 --- a/tests/pdegen/runs/amg_pde2d.inp +++ b/tests/pdegen/runs/amg_pde2d.inp @@ -46,10 +46,23 @@ FILTER ! Filtering of matrix: FILTER NOFILTER -2 ! Number of thresholds in vector, next line ignored if <= 0 0.05 0.025 ! Thresholds -0.0100d0 ! Smoothed aggregation threshold, ignored if < 0 -%%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% -BJAC ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC -ILU ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS SLU +%%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% +RKR ! Coarsest-level solver: MUMPS(global) UMF SLU SLUDIST JACOBI GS BJAC RKR(global) +NONE ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS(local) SLU RKR(local) +LLK ! AINV Variant for the coarse solver +BICGSTAB ! Krylov method for RKR solver/subsolver, ignored otherwise DIST ! Coarsest-level matrix distribution: DIST REPL -1 ! Coarsest-level fillin P for ILU(P) and ILU(T,P) +0 ! Coarsest-level fillin P for ILU(P) and ILU(T,P) +1 ! Coarsest-level inverse Fill level P for INVK 1.d-4 ! Coarsest-level threshold T for ILU(T,P) -1 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver +30 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver +30 ! maxit for RKR +1.d-4 ! eps for RKR +30 ! itrace for RKR +T ! Check the BJAC residual +T ! Print the BJAC residual +5 ! ITRACE for residual check +5 ! ITRACE for residual print +1.d-4 ! Tolerance for exit from BJAC +% dump +T ! Dump preconditioner diff --git a/tests/pdegen/runs/amg_pde3d.inp b/tests/pdegen/runs/amg_pde3d.inp index 26e3d719..82a1207f 100644 --- a/tests/pdegen/runs/amg_pde3d.inp +++ b/tests/pdegen/runs/amg_pde3d.inp @@ -46,10 +46,23 @@ NOFILTER ! Filtering of matrix: FILTER NOFILTER -2 ! Number of thresholds in vector, next line ignored if <= 0 0.05 0.025 ! Thresholds -0.0100d0 ! Smoothed aggregation threshold, ignored if < 0 -%%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% -BJAC ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC -ILU ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS SLU +%%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% +RKR ! Coarsest-level solver: MUMPS(global) UMF SLU SLUDIST JACOBI GS BJAC RKR(global) +NONE ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS(local) SLU RKR(local) +LLK ! AINV Variant for the coarse solver +BICGSTAB ! Krylov method for RKR solver/subsolver, ignored otherwise DIST ! Coarsest-level matrix distribution: DIST REPL -1 ! Coarsest-level fillin P for ILU(P) and ILU(T,P) +0 ! Coarsest-level fillin P for ILU(P) and ILU(T,P) +1 ! Coarsest-level inverse Fill level P for INVK 1.d-4 ! Coarsest-level threshold T for ILU(T,P) -1 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver +30 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver +30 ! maxit for RKR +1.d-4 ! eps for RKR +30 ! itrace for RKR +T ! Check the BJAC residual +T ! Print the BJAC residual +5 ! ITRACE for residual check +5 ! ITRACE for residual print +1.d-4 ! Tolerance for exit from BJAC +% dump +T ! Dump preconditioner