From 086d93dd28643177430890015a2414e2f7d30b22 Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Wed, 9 Dec 2020 15:45:02 +0100 Subject: [PATCH] RKR solver --- 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/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 ++++++++ 8 files changed, 3444 insertions(+) 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/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/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