From ea2f75776c7afabd345ccfc04f557476c44142ae Mon Sep 17 00:00:00 2001 From: sfilippone Date: Mon, 13 Nov 2023 16:58:30 +0100 Subject: [PATCH] Implement structure for polynomial smoother --- amgprec/Makefile | 2 + amgprec/amg_base_prec_type.F90 | 25 +- amgprec/amg_c_ilu_solver.f90 | 14 +- amgprec/amg_c_jac_solver.f90 | 5 +- amgprec/amg_d_beta_coeff_mod.f90 | 516 ++++++++++++++++++ amgprec/amg_d_ilu_solver.f90 | 14 +- amgprec/amg_d_jac_solver.f90 | 5 +- amgprec/amg_d_poly_smoother.f90 | 369 +++++++++++++ amgprec/amg_d_prec_mod.f90 | 1 + amgprec/amg_s_ilu_solver.f90 | 14 +- amgprec/amg_s_jac_solver.f90 | 5 +- amgprec/amg_z_ilu_solver.f90 | 14 +- amgprec/amg_z_jac_solver.f90 | 5 +- amgprec/impl/amg_c_smoothers_bld.f90 | 4 +- amgprec/impl/amg_ccprecset.F90 | 36 +- amgprec/impl/amg_d_smoothers_bld.f90 | 4 +- amgprec/impl/amg_dcprecset.F90 | 44 +- amgprec/impl/amg_dprecinit.F90 | 10 + amgprec/impl/amg_s_smoothers_bld.f90 | 4 +- amgprec/impl/amg_scprecset.F90 | 36 +- amgprec/impl/amg_z_smoothers_bld.f90 | 4 +- amgprec/impl/amg_zcprecset.F90 | 36 +- .../impl/level/amg_c_base_onelev_cseti.F90 | 2 +- .../impl/level/amg_d_base_onelev_csetc.F90 | 6 + .../impl/level/amg_d_base_onelev_cseti.F90 | 8 +- .../impl/level/amg_s_base_onelev_cseti.F90 | 2 +- .../impl/level/amg_z_base_onelev_cseti.F90 | 2 +- amgprec/impl/smoother/Makefile | 11 + .../smoother/amg_d_jac_smoother_clone.f90 | 1 + .../amg_d_poly_smoother_apply_vect.f90 | 357 ++++++++++++ .../impl/smoother/amg_d_poly_smoother_bld.f90 | 107 ++++ .../amg_d_poly_smoother_clear_data.f90 | 70 +++ .../smoother/amg_d_poly_smoother_clone.f90 | 90 +++ .../amg_d_poly_smoother_clone_settings.f90 | 96 ++++ .../impl/smoother/amg_d_poly_smoother_cnv.f90 | 77 +++ .../smoother/amg_d_poly_smoother_csetc.f90 | 72 +++ .../smoother/amg_d_poly_smoother_cseti.f90 | 69 +++ .../smoother/amg_d_poly_smoother_csetr.f90 | 69 +++ .../smoother/amg_d_poly_smoother_descr.f90 | 95 ++++ .../impl/smoother/amg_d_poly_smoother_dmp.f90 | 90 +++ amgprec/impl/solver/amg_c_ilu_solver_bld.f90 | 30 +- amgprec/impl/solver/amg_d_ilu_solver_bld.f90 | 30 +- amgprec/impl/solver/amg_s_ilu_solver_bld.f90 | 30 +- amgprec/impl/solver/amg_z_ilu_solver_bld.f90 | 30 +- samples/advanced/pdegen/amg_d_pde3d.F90 | 16 +- 45 files changed, 2364 insertions(+), 163 deletions(-) create mode 100644 amgprec/amg_d_beta_coeff_mod.f90 create mode 100644 amgprec/amg_d_poly_smoother.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_clone.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_cnv.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 create mode 100644 amgprec/impl/smoother/amg_d_poly_smoother_dmp.f90 diff --git a/amgprec/Makefile b/amgprec/Makefile index 81e6ce8c..0442f5ac 100644 --- a/amgprec/Makefile +++ b/amgprec/Makefile @@ -9,6 +9,7 @@ FINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(INCDIR) $(PSBLAS_INCLUDES) DMODOBJS=amg_d_prec_type.o \ amg_d_inner_mod.o amg_d_ilu_solver.o amg_d_diag_solver.o amg_d_jac_smoother.o amg_d_as_smoother.o \ + amg_d_poly_smoother.o amg_d_beta_coeff_mod.o\ amg_d_umf_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o amg_d_id_solver.o\ amg_d_base_solver_mod.o amg_d_base_smoother_mod.o amg_d_onelev_mod.o \ amg_d_gs_solver.o amg_d_mumps_solver.o amg_d_jac_solver.o \ @@ -164,6 +165,7 @@ amg_d_jac_smoother.o: amg_d_diag_solver.o amg_dprecinit.o amg_dprecset.o: amg_d_diag_solver.o amg_d_ilu_solver.o \ amg_d_umf_solver.o amg_d_as_smoother.o amg_d_jac_smoother.o \ amg_d_id_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o +amg_d_poly_smoother.o: amg_d_base_smoother_mod.o amg_d_beta_coeff_mod.o amg_s_mumps_solver.o amg_s_gs_solver.o amg_s_id_solver.o amg_s_slu_solver.o \ amg_s_diag_solver.o amg_s_ilu_solver.o amg_s_jac_solver.o: amg_s_base_solver_mod.o amg_s_prec_type.o diff --git a/amgprec/amg_base_prec_type.F90 b/amgprec/amg_base_prec_type.F90 index bd3ca19c..63f446c2 100644 --- a/amgprec/amg_base_prec_type.F90 +++ b/amgprec/amg_base_prec_type.F90 @@ -215,7 +215,8 @@ module amg_base_prec_type integer(psb_ipk_), parameter :: amg_fbgs_ = 6 integer(psb_ipk_), parameter :: amg_l1_gs_ = 7 integer(psb_ipk_), parameter :: amg_l1_fbgs_ = 8 - integer(psb_ipk_), parameter :: amg_max_prec_ = 8 + integer(psb_ipk_), parameter :: amg_poly_ = 9 + integer(psb_ipk_), parameter :: amg_max_prec_ = 9 ! ! Constants for pre/post signaling. Now only used internally ! @@ -233,9 +234,9 @@ module amg_base_prec_type integer(psb_ipk_), parameter :: amg_diag_scale_ = amg_slv_delta_+1 integer(psb_ipk_), parameter :: amg_l1_diag_scale_ = amg_slv_delta_+2 integer(psb_ipk_), parameter :: amg_gs_ = amg_slv_delta_+3 - ! !$ integer(psb_ipk_), parameter :: amg_ilu_n_ = amg_slv_delta_+4 - ! !$ integer(psb_ipk_), parameter :: amg_milu_n_ = amg_slv_delta_+5 - ! !$ integer(psb_ipk_), parameter :: amg_ilu_t_ = amg_slv_delta_+6 + integer(psb_ipk_), parameter :: amg_ilu_n_ = amg_slv_delta_+4 + integer(psb_ipk_), parameter :: amg_milu_n_ = amg_slv_delta_+5 + integer(psb_ipk_), parameter :: amg_ilu_t_ = amg_slv_delta_+6 integer(psb_ipk_), parameter :: amg_slu_ = amg_slv_delta_+7 integer(psb_ipk_), parameter :: amg_umf_ = amg_slv_delta_+8 integer(psb_ipk_), parameter :: amg_sludist_ = amg_slv_delta_+9 @@ -390,12 +391,12 @@ module amg_base_prec_type & ml_names(0:7)=(/'none ','additive ',& & 'multiplicative', 'VCycle ','WCycle ',& & 'KCycle ','KCycleSym ','new ML '/) - character(len=15), parameter :: & + character(len=16), parameter :: & & amg_fact_names(0:amg_max_sub_solve_)=(/& & 'none ','Jacobi ',& & 'L1-Jacobi ','none ','none ',& & 'none ','none ','L1-GS ',& - & 'L1-FBGS ','none ','Point Jacobi ',& + & 'L1-FBGS ','Polynomial ','none ','Point Jacobi ',& & 'L1-Jacobi ','Gauss-Seidel ','ILU(n) ',& & 'MILU(n) ','ILU(t,n) ',& & 'SuperLU ','UMFPACK LU ',& @@ -482,11 +483,11 @@ contains case('BGS','BWGS') val = amg_bwgs_ case('ILU') - val = psb_ilu_n_ + val = amg_ilu_n_ case('MILU') - val = psb_milu_n_ + val = amg_milu_n_ case('ILUT') - val = psb_ilu_t_ + val = amg_ilu_t_ case('MUMPS') val = amg_mumps_ case('UMF') @@ -557,6 +558,8 @@ contains val = amg_krm_ case('AS') val = amg_as_ + case('POLY') + val = amg_poly_ case('A_NORMI') val = amg_max_norm_ case('USER_CHOICE') @@ -1036,8 +1039,8 @@ contains integer(psb_ipk_), intent(in) :: ip logical :: is_legal_ilu_fact - is_legal_ilu_fact = ((ip==psb_ilu_n_).or.& - & (ip==psb_milu_n_).or.(ip==psb_ilu_t_)) + is_legal_ilu_fact = ((ip==amg_ilu_n_).or.& + & (ip==amg_milu_n_).or.(ip==amg_ilu_t_)) return end function is_legal_ilu_fact function is_legal_d_omega(ip) diff --git a/amgprec/amg_c_ilu_solver.f90 b/amgprec/amg_c_ilu_solver.f90 index 7a269d85..d9b4a1d2 100644 --- a/amgprec/amg_c_ilu_solver.f90 +++ b/amgprec/amg_c_ilu_solver.f90 @@ -234,7 +234,7 @@ contains ! Arguments class(amg_c_ilu_solver_type), intent(inout) :: sv - sv%fact_type = psb_ilu_n_ + sv%fact_type = amg_ilu_n_ sv%fill_in = 0 sv%thresh = szero @@ -255,13 +255,13 @@ contains info = psb_success_ call amg_check_def(sv%fact_type,& - & 'Factorization',psb_ilu_n_,is_legal_ilu_fact) + & 'Factorization',amg_ilu_n_,is_legal_ilu_fact) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) call amg_check_def(sv%fill_in,& & 'Level',izero,is_int_non_negative) - case(psb_ilu_t_) + case(amg_ilu_t_) call amg_check_def(sv%thresh,& & 'Eps',szero,is_legal_s_fact_thrs) end select @@ -439,9 +439,9 @@ contains write(iout_,*) trim(prefix_), ' Incomplete factorization solver: ',& & amg_fact_names(sv%fact_type) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in - case(psb_ilu_t_) + case(amg_ilu_t_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in write(iout_,*) trim(prefix_), ' Fill threshold :',sv%thresh end select @@ -496,7 +496,7 @@ contains implicit none integer(psb_ipk_) :: val - val = psb_ilu_n_ + val = amg_ilu_n_ end function c_ilu_solver_get_id function c_ilu_solver_get_wrksize() result(val) diff --git a/amgprec/amg_c_jac_solver.f90 b/amgprec/amg_c_jac_solver.f90 index 55335f1b..90ed85ce 100644 --- a/amgprec/amg_c_jac_solver.f90 +++ b/amgprec/amg_c_jac_solver.f90 @@ -403,7 +403,10 @@ contains info = psb_success_ call sv%a%free() - call sv%dv%free(info) + if (allocated(sv%dv)) then + call sv%dv%free(info) + deallocate(sv%dv) + end if if (allocated(sv%d)) deallocate(sv%d) call psb_erractionrestore(err_act) diff --git a/amgprec/amg_d_beta_coeff_mod.f90 b/amgprec/amg_d_beta_coeff_mod.f90 new file mode 100644 index 00000000..1bbeb876 --- /dev/null +++ b/amgprec/amg_d_beta_coeff_mod.f90 @@ -0,0 +1,516 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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_poly_smoother_mod.f90 +! +! Module: amg_d_poly_smoother_mod +! +! This module defines: +! the amg_d_poly_smoother_type data structure containing the +! smoother for a Jacobi/block Jacobi smoother. +! The smoother stores in ND the block off-diagonal matrix. +! One special case is treated separately, when the solver is DIAG or L1-DIAG +! then the ND is the entire off-diagonal part of the matrix (including the +! main diagonal block), so that it becomes possible to implement +! a pure Jacobi or L1-Jacobi global solver. +! +module amg_d_beta_coeff_mod + use psb_base_mod + + real(psb_dpk_), parameter :: amg_d_beta_vect(900) = [ & + & 1.1250000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0238728757031315_psb_dpk_, 1.2640890537108553_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0084254478202830_psb_dpk_, 1.0886783920873087_psb_dpk_, & + & 1.3375312590961856_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0039131042728535_psb_dpk_, 1.0403581118859304_psb_dpk_, & + & 1.1486349854625493_psb_dpk_, 1.3826886924100055_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0021293014616472_psb_dpk_, 1.0217371154926094_psb_dpk_, & + & 1.0787243319260302_psb_dpk_, 1.1981006529266300_psb_dpk_, & + & 1.4132254279168215_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0012851725594023_psb_dpk_, 1.0130429303523338_psb_dpk_, & + & 1.0467821512411335_psb_dpk_, 1.1161648941967548_psb_dpk_, & + & 1.2382902021844453_psb_dpk_, 1.4352429710674484_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0008346439791242_psb_dpk_, 1.0084394943012289_psb_dpk_, & + & 1.0300870776871385_psb_dpk_, 1.0740838409200377_psb_dpk_, & + & 1.1503618670736642_psb_dpk_, 1.2711647404613990_psb_dpk_, & + & 1.4518665864936395_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0005724663119766_psb_dpk_, 1.0057742766241562_psb_dpk_, & + & 1.0205018792294143_psb_dpk_, 1.0501980344456543_psb_dpk_, & + & 1.1011557298494106_psb_dpk_, 1.1808604280685657_psb_dpk_, & + & 1.2983858538257604_psb_dpk_, 1.4648607315109978_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0004096007283281_psb_dpk_, 1.0041243950610661_psb_dpk_, & + & 1.0146021214826659_psb_dpk_, 1.0356111362667175_psb_dpk_, & + & 1.0713997252919425_psb_dpk_, 1.1268827371096291_psb_dpk_, & + & 1.2078521914072933_psb_dpk_, 1.3212193071674674_psb_dpk_, & + & 1.4752964282069962_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0003031222965291_psb_dpk_, 1.0030484066079688_psb_dpk_, & + & 1.0107702271538761_psb_dpk_, 1.0261901159764004_psb_dpk_, & + & 1.0523172493375519_psb_dpk_, 1.0925574320754976_psb_dpk_, & + & 1.1508337666397197_psb_dpk_, 1.2317225087089441_psb_dpk_, & + & 1.3406080202445980_psb_dpk_, 1.4838612440701109_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0002305859520939_psb_dpk_, 1.0023167502402850_psb_dpk_, & + & 1.0081724539630488_psb_dpk_, 1.0198298656634219_psb_dpk_, & + & 1.0395021023532465_psb_dpk_, 1.0696504270054137_psb_dpk_, & + & 1.1130575429574259_psb_dpk_, 1.1729087627556418_psb_dpk_, & + & 1.2528830057679230_psb_dpk_, 1.3572557991951903_psb_dpk_, & + & 1.4910167256413891_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0001794720082837_psb_dpk_, 1.0018018913961957_psb_dpk_, & + & 1.0063486190730762_psb_dpk_, 1.0153786456630600_psb_dpk_, & + & 1.0305694283076039_psb_dpk_, 1.0537601969394355_psb_dpk_, & + & 1.0869986259207296_psb_dpk_, 1.1325918309791341_psb_dpk_, & + & 1.1931627335817252_psb_dpk_, 1.2717129367511055_psb_dpk_, & + & 1.3716933796979953_psb_dpk_, 1.4970841857556243_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0001424192155957_psb_dpk_, 1.0014290693262966_psb_dpk_, & + & 1.0050302898629815_psb_dpk_, 1.0121691051849540_psb_dpk_, & + & 1.0241487434279255_psb_dpk_, 1.0423815888082042_psb_dpk_, & + & 1.0684200812870084_psb_dpk_, 1.1039901093675994_psb_dpk_, & + & 1.1510274824264566_psb_dpk_, 1.2117181191012512_psb_dpk_, & + & 1.2885426486512805_psb_dpk_, 1.3843261938099158_psb_dpk_, & + & 1.5022941875736890_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0001149053826193_psb_dpk_, 1.0011524637691460_psb_dpk_, & + & 1.0040535733326481_psb_dpk_, 1.0097959057315313_psb_dpk_, & + & 1.0194130047299461_psb_dpk_, 1.0340142503543679_psb_dpk_, & + & 1.0548059960662932_psb_dpk_, 1.0831142030181304_psb_dpk_, & + & 1.1204089166089239_psb_dpk_, 1.1683309565544606_psb_dpk_, & + & 1.2287212228823874_psb_dpk_, 1.3036530570781755_psb_dpk_, & + & 1.3954681405367855_psb_dpk_, 1.5068164620958386_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000940475075257_psb_dpk_, 1.0009429169634352_psb_dpk_, & + & 1.0033144905644482_psb_dpk_, 1.0080029483381612_psb_dpk_, & + & 1.0158423625914039_psb_dpk_, 1.0277208331770495_psb_dpk_, & + & 1.0445953542283146_psb_dpk_, 1.0675076120612534_psb_dpk_, & + & 1.0976009254588965_psb_dpk_, 1.1361385536615733_psb_dpk_, & + & 1.1845236142623621_psb_dpk_, 1.2443208730447588_psb_dpk_, & + & 1.3172806908339272_psb_dpk_, 1.4053654389356023_psb_dpk_, & + & 1.5107787250184523_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000779482817921_psb_dpk_, 1.0007812684725339_psb_dpk_, & + & 1.0027448797440124_psb_dpk_, 1.0066229101701514_psb_dpk_, & + & 1.0130985883697137_psb_dpk_, 1.0228944832933697_psb_dpk_, & + & 1.0367832140998394_psb_dpk_, 1.0555987571989653_psb_dpk_, & + & 1.0802484840556024_psb_dpk_, 1.1117260713149764_psb_dpk_, & + & 1.1511254343107276_psb_dpk_, 1.1996558461497355_psb_dpk_, & + & 1.2586584174494597_psb_dpk_, 1.3296241265666493_psb_dpk_, & + & 1.4142136069557629_psb_dpk_, 1.5142789173034623_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000653242183546_psb_dpk_, 1.0006545722939437_psb_dpk_, & + & 1.0022987777448662_psb_dpk_, 1.0055432691173583_psb_dpk_, & + & 1.0109550075016893_psb_dpk_, 1.0191301541168694_psb_dpk_, & + & 1.0307019481191382_psb_dpk_, 1.0463489778000818_psb_dpk_, & + & 1.0668039321569163_psb_dpk_, 1.0928629244731740_psb_dpk_, & + & 1.1253954850882542_psb_dpk_, 1.1653553270075827_psb_dpk_, & + & 1.2137919954743157_psb_dpk_, 1.2718635211544003_psb_dpk_, & + & 1.3408502062615073_psb_dpk_, 1.4221696838526183_psb_dpk_, & + & 1.5173934027630227_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000552858792859_psb_dpk_, 1.0005538659610900_psb_dpk_, & + & 1.0019444166743086_psb_dpk_, 1.0046864301776393_psb_dpk_, & + & 1.0092557508630260_psb_dpk_, 1.0161502674772371_psb_dpk_, & + & 1.0258958148322650_psb_dpk_, 1.0390523408953256_psb_dpk_, & + & 1.0562203973533295_psb_dpk_, 1.0780480145522537_psb_dpk_, & + & 1.1052380250439366_psb_dpk_, 1.1385559038570177_psb_dpk_, & + & 1.1788381980793483_psb_dpk_, 1.2270016234308427_psb_dpk_, & + & 1.2840529112630572_psb_dpk_, 1.3510994958895055_psb_dpk_, & + & 1.4293611393851839_psb_dpk_, 1.5201825990516680_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000472036358790_psb_dpk_, 1.0004728102642675_psb_dpk_, & + & 1.0016593577469159_psb_dpk_, 1.0039976891368516_psb_dpk_, & + & 1.0078911941833455_psb_dpk_, 1.0137601583069535_psb_dpk_, & + & 1.0220462561721002_psb_dpk_, 1.0332172281153209_psb_dpk_, & + & 1.0477717791157513_psb_dpk_, 1.0662447417325256_psb_dpk_, & + & 1.0892125464929936_psb_dpk_, 1.1172990456131733_psb_dpk_, & + & 1.1511817386833911_psb_dpk_, 1.1915984520803475_psb_dpk_, & + & 1.2393545273929878_psb_dpk_, 1.2953305781018039_psb_dpk_, & + & 1.3604908781568688_psb_dpk_, 1.4358924509939206_psb_dpk_, & + & 1.5226949329440265_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000406232569254_psb_dpk_, 1.0004068351374691_psb_dpk_, & + & 1.0014274431564170_psb_dpk_, 1.0034377175807407_psb_dpk_, & + & 1.0067826854070978_psb_dpk_, 1.0118204999571436_psb_dpk_, & + & 1.0189259121271075_psb_dpk_, 1.0284938700470616_psb_dpk_, & + & 1.0409432748132981_psb_dpk_, 1.0567209210598594_psb_dpk_, & + & 1.0763056524407055_psb_dpk_, 1.1002127636100871_psb_dpk_, & + & 1.1289986820268283_psb_dpk_, 1.1632659648787138_psb_dpk_, & + & 1.2036686486408621_psb_dpk_, 1.2509179912601627_psb_dpk_, & + & 1.3057886497146727_psb_dpk_, 1.3691253387497200_psb_dpk_, & + & 1.4418500199624611_psb_dpk_, 1.5249696741164267_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000352114440929_psb_dpk_, 1.0003525892395289_psb_dpk_, & + & 1.0012368357172980_psb_dpk_, 1.0029777430511673_psb_dpk_, & + & 1.0058727830027672_psb_dpk_, 1.0102297507781717_psb_dpk_, & + & 1.0163694815733537_psb_dpk_, 1.0246286588536329_psb_dpk_, & + & 1.0353627340015590_psb_dpk_, 1.0489489776835172_psb_dpk_, & + & 1.0657896841306789_psb_dpk_, 1.0863155505114006_psb_dpk_, & + & 1.1109892546943501_psb_dpk_, 1.1403092559728156_psb_dpk_, & + & 1.1748138447471401_psb_dpk_, 1.2150854687543668_psb_dpk_, & + & 1.2617553651999671_psb_dpk_, 1.3155085300984379_psb_dpk_, & + & 1.3770890582780710_psb_dpk_, 1.4473058898645985_psb_dpk_, & + & 1.5270390016420912_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000307198714835_psb_dpk_, 1.0003075769178242_psb_dpk_, & + & 1.0010787281022711_psb_dpk_, 1.0025963829693492_psb_dpk_, & + & 1.0051188625231162_psb_dpk_, 1.0089126974249720_psb_dpk_, & + & 1.0142547789760521_psb_dpk_, 1.0214345766593154_psb_dpk_, & + & 1.0307564364069204_psb_dpk_, 1.0425419742322541_psb_dpk_, & + & 1.0571325804249445_psb_dpk_, 1.0748920501551993_psb_dpk_, & + & 1.0962093570737961_psb_dpk_, 1.1215015873309027_psb_dpk_, & + & 1.1512170523743910_psb_dpk_, 1.1858385999327761_psb_dpk_, & + & 1.2258871437439198_psb_dpk_, 1.2719254338660289_psb_dpk_, & + & 1.3245620908078453_psb_dpk_, 1.3844559282498121_psb_dpk_, & + & 1.4523205908039656_psb_dpk_, 1.5289295350887884_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000269609460124_psb_dpk_, 1.0002699137181752_psb_dpk_, & + & 1.0009464748475532_psb_dpk_, 1.0022775198638552_psb_dpk_, & + & 1.0044888368184179_psb_dpk_, 1.0078128087804721_psb_dpk_, & + & 1.0124901352066715_psb_dpk_, 1.0187716022931539_psb_dpk_, & + & 1.0269199126829005_psb_dpk_, 1.0372115852204526_psb_dpk_, & + & 1.0499389358225151_psb_dpk_, 1.0654121509688057_psb_dpk_, & + & 1.0839614658147161_psb_dpk_, 1.1059394594887115_psb_dpk_, & + & 1.1317234807654135_psb_dpk_, 1.1617182180038959_psb_dpk_, & + & 1.1963584280123116_psb_dpk_, 1.2361118393501820_psb_dpk_, & + & 1.2814822465106404_psb_dpk_, 1.3330128124440397_psb_dpk_, & + & 1.3912895979940381_psb_dpk_, 1.4569453380258381_psb_dpk_, & + & 1.5306634853375161_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000237911597230_psb_dpk_, 1.0002381585998457_psb_dpk_, & + & 1.0008349974382460_psb_dpk_, 1.0020088476285827_psb_dpk_, & + & 1.0039582343156432_psb_dpk_, 1.0068870298152559_psb_dpk_, & + & 1.0110058445931565_psb_dpk_, 1.0165334547611182_psb_dpk_, & + & 1.0236982737890488_psb_dpk_, 1.0327398763510158_psb_dpk_, & + & 1.0439105824804926_psb_dpk_, 1.0574771105088172_psb_dpk_, & + & 1.0737223076000839_psb_dpk_, 1.0929469670793606_psb_dpk_, & + & 1.1154717421787756_psb_dpk_, 1.1416391663018148_psb_dpk_, & + & 1.1718157904303341_psb_dpk_, 1.2063944488757254_psb_dpk_, & + & 1.2457966652063013_psb_dpk_, 1.2904752108716941_psb_dpk_, & + & 1.3409168297942540_psb_dpk_, 1.3976451430108305_psb_dpk_, & + & 1.4612237483301715_psb_dpk_, 1.5322595309246121_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000210994601235_psb_dpk_, 1.0002111968041199_psb_dpk_, & + & 1.0007403694573151_psb_dpk_, 1.0017808593384865_psb_dpk_, & + & 1.0035081686576977_psb_dpk_, 1.0061021720448531_psb_dpk_, & + & 1.0097482505685551_psb_dpk_, 1.0146384533048582_psb_dpk_, & + & 1.0209726922414943_psb_dpk_, 1.0289599764553270_psb_dpk_, & + & 1.0388196916802268_psb_dpk_, 1.0507829315895938_psb_dpk_, & + & 1.0650938873538003_psb_dpk_, 1.0820113022982043_psb_dpk_, & + & 1.1018099987843295_psb_dpk_, 1.1247824847650900_psb_dpk_, & + & 1.1512406478277994_psb_dpk_, 1.1815175449359154_psb_dpk_, & + & 1.2159692965153148_psb_dpk_, 1.2549770940040335_psb_dpk_, & + & 1.2989493304988182_psb_dpk_, 1.3483238646890843_psb_dpk_, & + & 1.4035704288718982_psb_dpk_, 1.4651931924923849_psb_dpk_, & + & 1.5337334933563860_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000187989989242_psb_dpk_, 1.0001881567984481_psb_dpk_, & + & 1.0006595227084085_psb_dpk_, 1.0015861311895899_psb_dpk_, & + & 1.0031239047778964_psb_dpk_, 1.0054323694760092_psb_dpk_, & + & 1.0086755868504005_psb_dpk_, 1.0130231071421940_psb_dpk_, & + & 1.0186509477893992_psb_dpk_, 1.0257426018654052_psb_dpk_, & + & 1.0344900810652515_psb_dpk_, 1.0450949980170887_psb_dpk_, & + & 1.0577696928624343_psb_dpk_, 1.0727384092356933_psb_dpk_, & + & 1.0902385249817814_psb_dpk_, 1.1105218431816117_psb_dpk_, & + & 1.1338559493090710_psb_dpk_, 1.1605256406217599_psb_dpk_, & + & 1.1908344341913664_psb_dpk_, 1.2251061603103259_psb_dpk_, & + & 1.2636866483695495_psb_dpk_, 1.3069455126904677_psb_dpk_, & + & 1.3552780462128098_psb_dpk_, 1.4091072303921326_psb_dpk_, & + & 1.4688858701459975_psb_dpk_, 1.5350988632115488_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000168211938973_psb_dpk_, 1.0001683505351420_psb_dpk_, & + & 1.0005900360142315_psb_dpk_, 1.0014188084960041_psb_dpk_, & + & 1.0027938311393803_psb_dpk_, 1.0048572584314193_psb_dpk_, & + & 1.0077550080990554_psb_dpk_, 1.0116375492127350_psb_dpk_, & + & 1.0166607098595459_psb_dpk_, 1.0229865078405374_psb_dpk_, & + & 1.0307840079371537_psb_dpk_, 1.0402302093961155_psb_dpk_, & + & 1.0515109674005423_psb_dpk_, 1.0648219524284319_psb_dpk_, & + & 1.0803696515480321_psb_dpk_, 1.0983724158638981_psb_dpk_, & + & 1.1190615585080472_psb_dpk_, 1.1426825077681895_psb_dpk_, & + & 1.1694960201606786_psb_dpk_, 1.1997794584895700_psb_dpk_, & + & 1.2338281401870808_psb_dpk_, 1.2719567615042522_psb_dpk_, & + & 1.3145009034164739_psb_dpk_, 1.3618186254259919_psb_dpk_, & + & 1.4142921537855777_psb_dpk_, 1.4723296710339275_psb_dpk_, & + & 1.5363672141264497_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000151113991291_psb_dpk_, 1.0001512299115287_psb_dpk_, & + & 1.0005299814085029_psb_dpk_, 1.0012742317597600_psb_dpk_, & + & 1.0025087130476142_psb_dpk_, 1.0043606572645858_psb_dpk_, & + & 1.0069604400315522_psb_dpk_, 1.0104422369100252_psb_dpk_, & + & 1.0149446949285030_psb_dpk_, 1.0206116219981500_psb_dpk_, & + & 1.0275926969588451_psb_dpk_, 1.0360442030716124_psb_dpk_, & + & 1.0461297878595799_psb_dpk_, 1.0580212522952626_psb_dpk_, & + & 1.0718993724396861_psb_dpk_, 1.0879547567564958_psb_dpk_, & + & 1.1063887424550545_psb_dpk_, 1.1274143343577541_psb_dpk_, & + & 1.1512571899424711_psb_dpk_, 1.1781566543781672_psb_dpk_, & + & 1.2083668495540898_psb_dpk_, 1.2421578212983135_psb_dpk_, & + & 1.2798167491932815_psb_dpk_, 1.3216492236219661_psb_dpk_, & + & 1.3679805949228399_psb_dpk_, 1.4191573997915068_psb_dpk_, & + & 1.4755488703473389_psb_dpk_, 1.5375485315807513_psb_dpk_, & + & 0.0000000000000000_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000136257096588_psb_dpk_, 1.0001363546506836_psb_dpk_, & + & 1.0004778107488095_psb_dpk_, 1.0011486612681773_psb_dpk_, & + & 1.0022611433613271_psb_dpk_, 1.0039295964948667_psb_dpk_, & + & 1.0062710027404669_psb_dpk_, 1.0094055369479136_psb_dpk_, & + & 1.0134571288503909_psb_dpk_, 1.0185540391932908_psb_dpk_, & + & 1.0248294520252528_psb_dpk_, 1.0324220853457433_psb_dpk_, & + & 1.0414768223656390_psb_dpk_, 1.0521453657079123_psb_dpk_, & + & 1.0645869169533493_psb_dpk_, 1.0789688840227822_psb_dpk_, & + & 1.0954676189818162_psb_dpk_, 1.1142691889576817_psb_dpk_, & + & 1.1355701829701565_psb_dpk_, 1.1595785576006521_psb_dpk_, & + & 1.1865145245551894_psb_dpk_, 1.2166114833191515_psb_dpk_, & + & 1.2501170022543431_psb_dpk_, 1.2872938516530203_psb_dpk_, & + & 1.3284210924391027_psb_dpk_, 1.3737952243949607_psb_dpk_, & + & 1.4237313979931023_psb_dpk_, 1.4785646941265451_psb_dpk_, & + & 1.5386514762605854_psb_dpk_, 0.0000000000000000_psb_dpk_, & + & 1.0000123285767939_psb_dpk_, 1.0001233683396147_psb_dpk_, & + & 1.0004322711781202_psb_dpk_, 1.0010390719329101_psb_dpk_, & + & 1.0020451337350940_psb_dpk_, 1.0035535979966428_psb_dpk_, & + & 1.0056698406248343_psb_dpk_, 1.0085019360540697_psb_dpk_, & + & 1.0121611307132341_psb_dpk_, 1.0167623275769953_psb_dpk_, & + & 1.0224245834847208_psb_dpk_, 1.0292716209515502_psb_dpk_, & + & 1.0374323562422998_psb_dpk_, 1.0470414455308106_psb_dpk_, & + & 1.0582398510249318_psb_dpk_, 1.0711754290010183_psb_dpk_, & + & 1.0860035417614331_psb_dpk_, 1.1028876956049132_psb_dpk_, & + & 1.1220002069820316_psb_dpk_, 1.1435228990979547_psb_dpk_, & + & 1.1676478313209715_psb_dpk_, 1.1945780638597872_psb_dpk_, & + & 1.2245284602839432_psb_dpk_, 1.2577265305821996_psb_dpk_, & + & 1.2944133175813315_psb_dpk_, 1.3348443296857557_psb_dpk_, & + & 1.3792905230439911_psb_dpk_, 1.4280393364047606_psb_dpk_, & + & 1.4813957820911738_psb_dpk_, 1.5396835966986973_psb_dpk_ ] + + + + +!!$ [1.1250000000000000_psb_dpk_, 0.0_psb_dpk_, 0.0_psb_dpk__psb_dpk_,,& +!!$ & 1.0238728757031315_psb_dpk_, 1.2640890537108553_psb_dpk_, 0.0_psb_dpk_,& +!!$ & 1.0084254478202830_psb_dpk_, 1.0886783920873087_psb_dpk_, 1.3375312590961856_psb_dpk_] + + real(psb_dpk_), parameter :: amg_d_beta_mat(30,30)=reshape(amg_d_beta_vect,[30,30]) + +end module amg_d_beta_coeff_mod diff --git a/amgprec/amg_d_ilu_solver.f90 b/amgprec/amg_d_ilu_solver.f90 index 00733655..c7813e6e 100644 --- a/amgprec/amg_d_ilu_solver.f90 +++ b/amgprec/amg_d_ilu_solver.f90 @@ -234,7 +234,7 @@ contains ! Arguments class(amg_d_ilu_solver_type), intent(inout) :: sv - sv%fact_type = psb_ilu_n_ + sv%fact_type = amg_ilu_n_ sv%fill_in = 0 sv%thresh = dzero @@ -255,13 +255,13 @@ contains info = psb_success_ call amg_check_def(sv%fact_type,& - & 'Factorization',psb_ilu_n_,is_legal_ilu_fact) + & 'Factorization',amg_ilu_n_,is_legal_ilu_fact) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) call amg_check_def(sv%fill_in,& & 'Level',izero,is_int_non_negative) - case(psb_ilu_t_) + case(amg_ilu_t_) call amg_check_def(sv%thresh,& & 'Eps',dzero,is_legal_d_fact_thrs) end select @@ -439,9 +439,9 @@ contains write(iout_,*) trim(prefix_), ' Incomplete factorization solver: ',& & amg_fact_names(sv%fact_type) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in - case(psb_ilu_t_) + case(amg_ilu_t_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in write(iout_,*) trim(prefix_), ' Fill threshold :',sv%thresh end select @@ -496,7 +496,7 @@ contains implicit none integer(psb_ipk_) :: val - val = psb_ilu_n_ + val = amg_ilu_n_ end function d_ilu_solver_get_id function d_ilu_solver_get_wrksize() result(val) diff --git a/amgprec/amg_d_jac_solver.f90 b/amgprec/amg_d_jac_solver.f90 index 25bb1375..eb7c93ce 100644 --- a/amgprec/amg_d_jac_solver.f90 +++ b/amgprec/amg_d_jac_solver.f90 @@ -403,7 +403,10 @@ contains info = psb_success_ call sv%a%free() - call sv%dv%free(info) + if (allocated(sv%dv)) then + call sv%dv%free(info) + deallocate(sv%dv) + end if if (allocated(sv%d)) deallocate(sv%d) call psb_erractionrestore(err_act) diff --git a/amgprec/amg_d_poly_smoother.f90 b/amgprec/amg_d_poly_smoother.f90 new file mode 100644 index 00000000..aa1dd5d1 --- /dev/null +++ b/amgprec/amg_d_poly_smoother.f90 @@ -0,0 +1,369 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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_poly_smoother_mod.f90 +! +! Module: amg_d_poly_smoother_mod +! +! This module defines: +! the amg_d_poly_smoother_type data structure containing the +! smoother for a Jacobi/block Jacobi smoother. +! The smoother stores in ND the block off-diagonal matrix. +! One special case is treated separately, when the solver is DIAG or L1-DIAG +! then the ND is the entire off-diagonal part of the matrix (including the +! main diagonal block), so that it becomes possible to implement +! a pure Jacobi or L1-Jacobi global solver. +! +module amg_d_poly_smoother + use amg_d_base_smoother_mod + use amg_d_beta_coeff_mod + + type, extends(amg_d_base_smoother_type) :: amg_d_poly_smoother_type + ! The local solver component is inherited from the + ! parent type. + ! class(amg_d_base_solver_type), allocatable :: sv + ! + integer(psb_ipk_) :: pdegree + type(psb_dspmat_type), pointer :: pa => null() + real(psb_dpk_), allocatable :: poly_beta(:) + real(psb_dpk_) :: rho_ba + contains + procedure, pass(sm) :: apply_v => amg_d_poly_smoother_apply_vect +!!$ procedure, pass(sm) :: apply_a => amg_d_poly_smoother_apply + procedure, pass(sm) :: dump => amg_d_poly_smoother_dmp + procedure, pass(sm) :: build => amg_d_poly_smoother_bld + procedure, pass(sm) :: cnv => amg_d_poly_smoother_cnv + procedure, pass(sm) :: clone => amg_d_poly_smoother_clone + procedure, pass(sm) :: clone_settings => amg_d_poly_smoother_clone_settings + procedure, pass(sm) :: clear_data => amg_d_poly_smoother_clear_data + procedure, pass(sm) :: free => d_poly_smoother_free + procedure, pass(sm) :: cseti => amg_d_poly_smoother_cseti + procedure, pass(sm) :: csetc => amg_d_poly_smoother_csetc + procedure, pass(sm) :: csetr => amg_d_poly_smoother_csetr + procedure, pass(sm) :: descr => amg_d_poly_smoother_descr + procedure, pass(sm) :: sizeof => d_poly_smoother_sizeof + procedure, pass(sm) :: default => d_poly_smoother_default + procedure, pass(sm) :: get_nzeros => d_poly_smoother_get_nzeros + procedure, pass(sm) :: get_wrksz => d_poly_smoother_get_wrksize + procedure, nopass :: get_fmt => d_poly_smoother_get_fmt + procedure, nopass :: get_id => d_poly_smoother_get_id + end type amg_d_poly_smoother_type + private :: d_poly_smoother_free, & + & d_poly_smoother_sizeof, d_poly_smoother_get_nzeros, & + & d_poly_smoother_get_fmt, d_poly_smoother_get_id, & + & d_poly_smoother_get_wrksize + + + interface + subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& + & sweeps,work,wv,info,init,initu) + import :: psb_desc_type, amg_d_poly_smoother_type, psb_d_vect_type, psb_dpk_, & + & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,& + & psb_ipk_ + + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_poly_smoother_type), intent(inout) :: sm + 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 + integer(psb_ipk_), intent(in) :: sweeps + 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_poly_smoother_apply_vect + end interface + +!!$ interface +!!$ subroutine amg_d_poly_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& +!!$ & sweeps,work,info,init,initu) +!!$ import :: psb_desc_type, amg_d_poly_smoother_type, psb_d_vect_type, psb_dpk_, & +!!$ & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, & +!!$ & psb_ipk_ +!!$ type(psb_desc_type), intent(in) :: desc_data +!!$ class(amg_d_poly_smoother_type), intent(inout) :: sm +!!$ 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 +!!$ integer(psb_ipk_), intent(in) :: sweeps +!!$ 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_poly_smoother_apply +!!$ end interface +!!$ + + interface + subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) + import :: psb_desc_type, amg_d_poly_smoother_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 + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_d_poly_smoother_bld + end interface + + interface + subroutine amg_d_poly_smoother_cnv(sm,info,amold,vmold,imold) + import :: amg_d_poly_smoother_type, psb_dpk_, & + & psb_d_base_sparse_mat, psb_d_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + class(amg_d_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine amg_d_poly_smoother_cnv + end interface + + interface + subroutine amg_d_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solver,global_num) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, amg_d_poly_smoother_type, psb_epk_, psb_desc_type, & + & psb_ipk_ + implicit none + class(amg_d_poly_smoother_type), intent(in) :: sm + 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) :: smoother, solver, global_num + end subroutine amg_d_poly_smoother_dmp + end interface + + interface + subroutine amg_d_poly_smoother_clone(sm,smout,info) + import :: amg_d_poly_smoother_type, psb_dpk_, & + & amg_d_base_smoother_type, psb_ipk_ + class(amg_d_poly_smoother_type), intent(inout) :: sm + class(amg_d_base_smoother_type), allocatable, intent(inout) :: smout + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_poly_smoother_clone + end interface + + interface + subroutine amg_d_poly_smoother_clone_settings(sm,smout,info) + import :: amg_d_poly_smoother_type, psb_dpk_, & + & amg_d_base_smoother_type, psb_ipk_ + class(amg_d_poly_smoother_type), intent(inout) :: sm + class(amg_d_base_smoother_type), allocatable, intent(inout) :: smout + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_poly_smoother_clone_settings + end interface + + interface + subroutine amg_d_poly_smoother_clear_data(sm,info) + import :: amg_d_poly_smoother_type, psb_dpk_, & + & amg_d_base_smoother_type, psb_ipk_ + class(amg_d_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_poly_smoother_clear_data + end interface + + interface + subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) + import :: amg_d_poly_smoother_type, psb_ipk_ + class(amg_d_poly_smoother_type), intent(in) :: sm + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + character(len=*), intent(in), optional :: prefix + end subroutine amg_d_poly_smoother_descr + end interface + + interface + subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, amg_d_poly_smoother_type, psb_epk_, psb_desc_type, psb_ipk_ + implicit none + class(amg_d_poly_smoother_type), intent(inout) :: sm + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + end subroutine amg_d_poly_smoother_cseti + end interface + + interface + subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, amg_d_poly_smoother_type, psb_epk_, psb_desc_type, psb_ipk_ + implicit none + class(amg_d_poly_smoother_type), intent(inout) :: sm + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + end subroutine amg_d_poly_smoother_csetc + end interface + + interface + subroutine amg_d_poly_smoother_csetr(sm,what,val,info,idx) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, amg_d_poly_smoother_type, psb_epk_, psb_desc_type, psb_ipk_ + implicit none + class(amg_d_poly_smoother_type), intent(inout) :: sm + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + end subroutine amg_d_poly_smoother_csetr + end interface + + +contains + + + subroutine d_poly_smoother_free(sm,info) + + + Implicit None + + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_poly_smoother_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + + + if (allocated(sm%sv)) then + call sm%sv%free(info) + if (info == psb_success_) deallocate(sm%sv,stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + end if + if (allocated(sm%poly_beta)) deallocate(sm%poly_beta) + sm%pa => null() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_poly_smoother_free + + function d_poly_smoother_sizeof(sm) result(val) + + implicit none + ! Arguments + class(amg_d_poly_smoother_type), intent(in) :: sm + integer(psb_epk_) :: val + + val = psb_sizeof_dp + if (allocated(sm%sv)) val = val + sm%sv%sizeof() + if (allocated(sm%poly_beta)) val = val + psb_sizeof_dp * size(sm%poly_beta) + + return + end function d_poly_smoother_sizeof + + subroutine d_poly_smoother_default(sm) + + Implicit None + + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + + ! + ! Default: BJAC with no residual check + ! + sm%pdegree = 1 + sm%rho_ba = dzero + + if (allocated(sm%sv)) then + call sm%sv%default() + end if + + return + end subroutine d_poly_smoother_default + + function d_poly_smoother_get_nzeros(sm) result(val) + + implicit none + ! Arguments + class(amg_d_poly_smoother_type), intent(in) :: sm + integer(psb_epk_) :: val + integer(psb_ipk_) :: i + + val = 0 + if (allocated(sm%sv)) val = val + sm%sv%get_nzeros() + + return + end function d_poly_smoother_get_nzeros + + function d_poly_smoother_get_wrksize(sm) result(val) + implicit none + class(amg_d_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_) :: val + + val = 4 + if (allocated(sm%sv)) val = val + sm%sv%get_wrksz() + + end function d_poly_smoother_get_wrksize + + function d_poly_smoother_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Polynomial smoother" + end function d_poly_smoother_get_fmt + + function d_poly_smoother_get_id() result(val) + implicit none + integer(psb_ipk_) :: val + + val = amg_poly_ + end function d_poly_smoother_get_id + + +end module amg_d_poly_smoother diff --git a/amgprec/amg_d_prec_mod.f90 b/amgprec/amg_d_prec_mod.f90 index 96e2f974..f86133d8 100644 --- a/amgprec/amg_d_prec_mod.f90 +++ b/amgprec/amg_d_prec_mod.f90 @@ -47,6 +47,7 @@ module amg_d_prec_mod use amg_d_prec_type use amg_d_jac_smoother use amg_d_as_smoother + use amg_d_poly_smoother use amg_d_id_solver use amg_d_diag_solver use amg_d_l1_diag_solver diff --git a/amgprec/amg_s_ilu_solver.f90 b/amgprec/amg_s_ilu_solver.f90 index dd642746..8bca532b 100644 --- a/amgprec/amg_s_ilu_solver.f90 +++ b/amgprec/amg_s_ilu_solver.f90 @@ -234,7 +234,7 @@ contains ! Arguments class(amg_s_ilu_solver_type), intent(inout) :: sv - sv%fact_type = psb_ilu_n_ + sv%fact_type = amg_ilu_n_ sv%fill_in = 0 sv%thresh = szero @@ -255,13 +255,13 @@ contains info = psb_success_ call amg_check_def(sv%fact_type,& - & 'Factorization',psb_ilu_n_,is_legal_ilu_fact) + & 'Factorization',amg_ilu_n_,is_legal_ilu_fact) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) call amg_check_def(sv%fill_in,& & 'Level',izero,is_int_non_negative) - case(psb_ilu_t_) + case(amg_ilu_t_) call amg_check_def(sv%thresh,& & 'Eps',szero,is_legal_s_fact_thrs) end select @@ -439,9 +439,9 @@ contains write(iout_,*) trim(prefix_), ' Incomplete factorization solver: ',& & amg_fact_names(sv%fact_type) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in - case(psb_ilu_t_) + case(amg_ilu_t_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in write(iout_,*) trim(prefix_), ' Fill threshold :',sv%thresh end select @@ -496,7 +496,7 @@ contains implicit none integer(psb_ipk_) :: val - val = psb_ilu_n_ + val = amg_ilu_n_ end function s_ilu_solver_get_id function s_ilu_solver_get_wrksize() result(val) diff --git a/amgprec/amg_s_jac_solver.f90 b/amgprec/amg_s_jac_solver.f90 index 8cc66bbc..0ecbd10d 100644 --- a/amgprec/amg_s_jac_solver.f90 +++ b/amgprec/amg_s_jac_solver.f90 @@ -403,7 +403,10 @@ contains info = psb_success_ call sv%a%free() - call sv%dv%free(info) + if (allocated(sv%dv)) then + call sv%dv%free(info) + deallocate(sv%dv) + end if if (allocated(sv%d)) deallocate(sv%d) call psb_erractionrestore(err_act) diff --git a/amgprec/amg_z_ilu_solver.f90 b/amgprec/amg_z_ilu_solver.f90 index 48b5ff1f..f2c6f29b 100644 --- a/amgprec/amg_z_ilu_solver.f90 +++ b/amgprec/amg_z_ilu_solver.f90 @@ -234,7 +234,7 @@ contains ! Arguments class(amg_z_ilu_solver_type), intent(inout) :: sv - sv%fact_type = psb_ilu_n_ + sv%fact_type = amg_ilu_n_ sv%fill_in = 0 sv%thresh = dzero @@ -255,13 +255,13 @@ contains info = psb_success_ call amg_check_def(sv%fact_type,& - & 'Factorization',psb_ilu_n_,is_legal_ilu_fact) + & 'Factorization',amg_ilu_n_,is_legal_ilu_fact) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) call amg_check_def(sv%fill_in,& & 'Level',izero,is_int_non_negative) - case(psb_ilu_t_) + case(amg_ilu_t_) call amg_check_def(sv%thresh,& & 'Eps',dzero,is_legal_d_fact_thrs) end select @@ -439,9 +439,9 @@ contains write(iout_,*) trim(prefix_), ' Incomplete factorization solver: ',& & amg_fact_names(sv%fact_type) select case(sv%fact_type) - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in - case(psb_ilu_t_) + case(amg_ilu_t_) write(iout_,*) trim(prefix_), ' Fill level:',sv%fill_in write(iout_,*) trim(prefix_), ' Fill threshold :',sv%thresh end select @@ -496,7 +496,7 @@ contains implicit none integer(psb_ipk_) :: val - val = psb_ilu_n_ + val = amg_ilu_n_ end function z_ilu_solver_get_id function z_ilu_solver_get_wrksize() result(val) diff --git a/amgprec/amg_z_jac_solver.f90 b/amgprec/amg_z_jac_solver.f90 index 28f1199e..5d273537 100644 --- a/amgprec/amg_z_jac_solver.f90 +++ b/amgprec/amg_z_jac_solver.f90 @@ -403,7 +403,10 @@ contains info = psb_success_ call sv%a%free() - call sv%dv%free(info) + if (allocated(sv%dv)) then + call sv%dv%free(info) + deallocate(sv%dv) + end if if (allocated(sv%d)) deallocate(sv%d) call psb_erractionrestore(err_act) diff --git a/amgprec/impl/amg_c_smoothers_bld.f90 b/amgprec/impl/amg_c_smoothers_bld.f90 index 8ad4d6eb..d684e389 100644 --- a/amgprec/impl/amg_c_smoothers_bld.f90 +++ b/amgprec/impl/amg_c_smoothers_bld.f90 @@ -186,8 +186,8 @@ subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) & ' but it has been changed to distributed.' end if - case(psb_ilu_n_, psb_ilu_t_,psb_milu_n_) - if (prec%precv(iszv)%sm%sv%get_id() /= psb_ilu_n_) then + case(amg_ilu_n_, amg_ilu_t_,amg_milu_n_) + if (prec%precv(iszv)%sm%sv%get_id() /= amg_ilu_n_) then write(psb_err_unit,*) & & 'AMG4PSBLAS: Warning: original coarse solver was requested as ',& & amg_fact_names(coarse_solve_id) diff --git a/amgprec/impl/amg_ccprecset.F90 b/amgprec/impl/amg_ccprecset.F90 index 70dc3013..4e7c7c4e 100644 --- a/amgprec/impl/amg_ccprecset.F90 +++ b/amgprec/impl/amg_ccprecset.F90 @@ -471,7 +471,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -485,7 +485,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -508,7 +508,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -516,21 +516,21 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -545,7 +545,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -568,7 +568,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -592,7 +592,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -700,7 +700,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -713,7 +713,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -736,7 +736,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -744,21 +744,21 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -773,7 +773,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -796,7 +796,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -820,7 +820,7 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) diff --git a/amgprec/impl/amg_d_smoothers_bld.f90 b/amgprec/impl/amg_d_smoothers_bld.f90 index 76347dc4..c2bf3b99 100644 --- a/amgprec/impl/amg_d_smoothers_bld.f90 +++ b/amgprec/impl/amg_d_smoothers_bld.f90 @@ -186,8 +186,8 @@ subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) & ' but it has been changed to distributed.' end if - case(psb_ilu_n_, psb_ilu_t_,psb_milu_n_) - if (prec%precv(iszv)%sm%sv%get_id() /= psb_ilu_n_) then + case(amg_ilu_n_, amg_ilu_t_,amg_milu_n_) + if (prec%precv(iszv)%sm%sv%get_id() /= amg_ilu_n_) then write(psb_err_unit,*) & & 'AMG4PSBLAS: Warning: original coarse solver was requested as ',& & amg_fact_names(coarse_solve_id) diff --git a/amgprec/impl/amg_dcprecset.F90 b/amgprec/impl/amg_dcprecset.F90 index deaced0d..c46f5224 100644 --- a/amgprec/impl/amg_dcprecset.F90 +++ b/amgprec/impl/amg_dcprecset.F90 @@ -81,6 +81,7 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx) use amg_d_prec_mod, amg_protect_name => amg_dcprecseti use amg_d_jac_smoother use amg_d_as_smoother + use amg_d_poly_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_ilu_solver @@ -125,7 +126,7 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx) info = 3111 write(psb_err_unit,*) name,& & ': Error: uninitialized preconditioner,',& - &' should call amg_PRECINIT' + & ' should call amg_PRECINIT' return endif @@ -312,6 +313,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) use amg_d_prec_mod, amg_protect_name => amg_dcprecsetc use amg_d_jac_smoother use amg_d_as_smoother + use amg_d_poly_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_ilu_solver @@ -402,6 +404,10 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) do il=ilev_, ilmax_ call p%precv(il)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) end do + case ('POLY') + do il=ilev_, ilmax_ + call p%precv(il)%set('SMOOTHER_TYPE',amg_poly_,info,pos=pos) + end do case ('L1-BJAC') do il=ilev_, ilmax_ call p%precv(il)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos) @@ -485,7 +491,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -501,7 +507,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -524,7 +530,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -532,21 +538,21 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -561,7 +567,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -591,7 +597,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -629,7 +635,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -739,7 +745,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -754,7 +760,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -777,7 +783,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -785,21 +791,21 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -814,7 +820,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -844,7 +850,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -882,7 +888,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) diff --git a/amgprec/impl/amg_dprecinit.F90 b/amgprec/impl/amg_dprecinit.F90 index eaa861e9..491d8eca 100644 --- a/amgprec/impl/amg_dprecinit.F90 +++ b/amgprec/impl/amg_dprecinit.F90 @@ -93,6 +93,7 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) use psb_base_mod use amg_d_prec_mod, amg_protect_name => amg_dprecinit use amg_d_jac_smoother + use amg_d_poly_smoother use amg_d_as_smoother use amg_d_id_solver use amg_d_diag_solver @@ -156,6 +157,15 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() + case ('POLY') + nlev_ = 1 + ilev_ = 1 + allocate(prec%precv(nlev_),stat=info) + allocate(amg_d_poly_smoother_type :: prec%precv(ilev_)%sm, stat=info) + if (info /= psb_success_) return + allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) + call prec%precv(ilev_)%default() + case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') nlev_ = 1 ilev_ = 1 diff --git a/amgprec/impl/amg_s_smoothers_bld.f90 b/amgprec/impl/amg_s_smoothers_bld.f90 index 8149d3bb..00385c8a 100644 --- a/amgprec/impl/amg_s_smoothers_bld.f90 +++ b/amgprec/impl/amg_s_smoothers_bld.f90 @@ -186,8 +186,8 @@ subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) & ' but it has been changed to distributed.' end if - case(psb_ilu_n_, psb_ilu_t_,psb_milu_n_) - if (prec%precv(iszv)%sm%sv%get_id() /= psb_ilu_n_) then + case(amg_ilu_n_, amg_ilu_t_,amg_milu_n_) + if (prec%precv(iszv)%sm%sv%get_id() /= amg_ilu_n_) then write(psb_err_unit,*) & & 'AMG4PSBLAS: Warning: original coarse solver was requested as ',& & amg_fact_names(coarse_solve_id) diff --git a/amgprec/impl/amg_scprecset.F90 b/amgprec/impl/amg_scprecset.F90 index 754725ad..741325fb 100644 --- a/amgprec/impl/amg_scprecset.F90 +++ b/amgprec/impl/amg_scprecset.F90 @@ -471,7 +471,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -485,7 +485,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -508,7 +508,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -516,21 +516,21 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -545,7 +545,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -568,7 +568,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -592,7 +592,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -700,7 +700,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -713,7 +713,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -736,7 +736,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -744,21 +744,21 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -773,7 +773,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -796,7 +796,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -820,7 +820,7 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) diff --git a/amgprec/impl/amg_z_smoothers_bld.f90 b/amgprec/impl/amg_z_smoothers_bld.f90 index 95293993..eeaef05b 100644 --- a/amgprec/impl/amg_z_smoothers_bld.f90 +++ b/amgprec/impl/amg_z_smoothers_bld.f90 @@ -186,8 +186,8 @@ subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) & ' but it has been changed to distributed.' end if - case(psb_ilu_n_, psb_ilu_t_,psb_milu_n_) - if (prec%precv(iszv)%sm%sv%get_id() /= psb_ilu_n_) then + case(amg_ilu_n_, amg_ilu_t_,amg_milu_n_) + if (prec%precv(iszv)%sm%sv%get_id() /= amg_ilu_n_) then write(psb_err_unit,*) & & 'AMG4PSBLAS: Warning: original coarse solver was requested as ',& & amg_fact_names(coarse_solve_id) diff --git a/amgprec/impl/amg_zcprecset.F90 b/amgprec/impl/amg_zcprecset.F90 index 5317ca6c..3b63befd 100644 --- a/amgprec/impl/amg_zcprecset.F90 +++ b/amgprec/impl/amg_zcprecset.F90 @@ -485,7 +485,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -501,7 +501,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -524,7 +524,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -532,21 +532,21 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -561,7 +561,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -591,7 +591,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -629,7 +629,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -739,7 +739,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -754,7 +754,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #elif defined(HAVE_MUMPS_) call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos) #else - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) #endif if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& @@ -777,7 +777,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -785,21 +785,21 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) #endif case('ILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('ILUT') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_t_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_t_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) call p%precv(nlev_)%set('COARSE_MAT',amg_repl_mat_,info,pos=pos) case('MILU') call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_milu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_milu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_repl_mat_) @@ -814,7 +814,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -844,7 +844,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) @@ -882,7 +882,7 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info,pos=pos) #else call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) - call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos) + call p%precv(nlev_)%set('SUB_SOLVE',amg_ilu_n_,info,pos=pos) if (hier_asb) & & call amg_warn_coarse_mat(p%precv(nlev_)%parms%get_coarse_mat(),& & amg_distr_mat_) diff --git a/amgprec/impl/level/amg_c_base_onelev_cseti.F90 b/amgprec/impl/level/amg_c_base_onelev_cseti.F90 index deba9001..1cc22cd9 100644 --- a/amgprec/impl/level/amg_c_base_onelev_cseti.F90 +++ b/amgprec/impl/level/amg_c_base_onelev_cseti.F90 @@ -164,7 +164,7 @@ subroutine amg_c_base_onelev_cseti(lv,what,val,info,pos,idx) case (amg_bwgs_) call lv%set(amg_c_bwgs_solver_mold,info,pos=pos) - case (psb_ilu_n_,psb_milu_n_,psb_ilu_t_) + case (amg_ilu_n_,amg_milu_n_,amg_ilu_t_) call lv%set(amg_c_ilu_solver_mold,info,pos=pos) if (info == 0) then if ((ipos_==amg_smooth_pre_) .or.(ipos_==amg_smooth_both_)) then diff --git a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 index 94f45276..c8a8941d 100644 --- a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 @@ -45,6 +45,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) use amg_d_parmatch_aggregator_mod use amg_d_jac_smoother use amg_d_as_smoother + use amg_d_poly_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_jac_solver @@ -84,6 +85,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) type(amg_d_jac_smoother_type) :: amg_d_jac_smoother_mold type(amg_d_l1_jac_smoother_type) :: amg_d_l1_jac_smoother_mold type(amg_d_as_smoother_type) :: amg_d_as_smoother_mold + type(amg_d_poly_smoother_type) :: amg_d_poly_smoother_mold type(amg_d_diag_solver_type) :: amg_d_diag_solver_mold type(amg_d_l1_diag_solver_type) :: amg_d_l1_diag_solver_mold type(amg_d_jac_solver_type) :: amg_d_jac_solver_mold @@ -156,6 +158,10 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) call lv%set(amg_d_as_smoother_mold,info,pos=pos) if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos) + case ('POLY') + call lv%set(amg_d_poly_smoother_mold,info,pos=pos) + if (info == 0) call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos) + case ('GS','FWGS') call lv%set(amg_d_jac_smoother_mold,info,pos='pre') if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre') diff --git a/amgprec/impl/level/amg_d_base_onelev_cseti.F90 b/amgprec/impl/level/amg_d_base_onelev_cseti.F90 index b5ca549b..195fbed7 100644 --- a/amgprec/impl/level/amg_d_base_onelev_cseti.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_cseti.F90 @@ -45,6 +45,7 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx) use amg_d_parmatch_aggregator_mod use amg_d_jac_smoother use amg_d_as_smoother + use amg_d_poly_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_ilu_solver @@ -79,6 +80,7 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx) type(amg_d_jac_smoother_type) :: amg_d_jac_smoother_mold type(amg_d_l1_jac_smoother_type) :: amg_d_l1_jac_smoother_mold type(amg_d_as_smoother_type) :: amg_d_as_smoother_mold + type(amg_d_poly_smoother_type) :: amg_d_poly_smoother_mold type(amg_d_diag_solver_type) :: amg_d_diag_solver_mold type(amg_d_l1_diag_solver_type) :: amg_d_l1_diag_solver_mold type(amg_d_ilu_solver_type) :: amg_d_ilu_solver_mold @@ -141,6 +143,10 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx) call lv%set(amg_d_as_smoother_mold,info,pos=pos) if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos) + case (amg_poly_) + call lv%set(amg_d_poly_smoother_mold,info,pos=pos) + if (info == 0) call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos) + case (amg_fbgs_) call lv%set(amg_d_jac_smoother_mold,info,pos='pre') if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre') @@ -177,7 +183,7 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx) case (amg_bwgs_) call lv%set(amg_d_bwgs_solver_mold,info,pos=pos) - case (psb_ilu_n_,psb_milu_n_,psb_ilu_t_) + case (amg_ilu_n_,amg_milu_n_,amg_ilu_t_) call lv%set(amg_d_ilu_solver_mold,info,pos=pos) if (info == 0) then if ((ipos_==amg_smooth_pre_) .or.(ipos_==amg_smooth_both_)) then diff --git a/amgprec/impl/level/amg_s_base_onelev_cseti.F90 b/amgprec/impl/level/amg_s_base_onelev_cseti.F90 index 1211a662..74fb3899 100644 --- a/amgprec/impl/level/amg_s_base_onelev_cseti.F90 +++ b/amgprec/impl/level/amg_s_base_onelev_cseti.F90 @@ -165,7 +165,7 @@ subroutine amg_s_base_onelev_cseti(lv,what,val,info,pos,idx) case (amg_bwgs_) call lv%set(amg_s_bwgs_solver_mold,info,pos=pos) - case (psb_ilu_n_,psb_milu_n_,psb_ilu_t_) + case (amg_ilu_n_,amg_milu_n_,amg_ilu_t_) call lv%set(amg_s_ilu_solver_mold,info,pos=pos) if (info == 0) then if ((ipos_==amg_smooth_pre_) .or.(ipos_==amg_smooth_both_)) then diff --git a/amgprec/impl/level/amg_z_base_onelev_cseti.F90 b/amgprec/impl/level/amg_z_base_onelev_cseti.F90 index b6a447a4..f50a785f 100644 --- a/amgprec/impl/level/amg_z_base_onelev_cseti.F90 +++ b/amgprec/impl/level/amg_z_base_onelev_cseti.F90 @@ -176,7 +176,7 @@ subroutine amg_z_base_onelev_cseti(lv,what,val,info,pos,idx) case (amg_bwgs_) call lv%set(amg_z_bwgs_solver_mold,info,pos=pos) - case (psb_ilu_n_,psb_milu_n_,psb_ilu_t_) + case (amg_ilu_n_,amg_milu_n_,amg_ilu_t_) call lv%set(amg_z_ilu_solver_mold,info,pos=pos) if (info == 0) then if ((ipos_==amg_smooth_pre_) .or.(ipos_==amg_smooth_both_)) then diff --git a/amgprec/impl/smoother/Makefile b/amgprec/impl/smoother/Makefile index f26b8f00..58884bc8 100644 --- a/amgprec/impl/smoother/Makefile +++ b/amgprec/impl/smoother/Makefile @@ -97,6 +97,17 @@ amg_d_jac_smoother_csetr.o \ amg_d_l1_jac_smoother_bld.o \ amg_d_l1_jac_smoother_descr.o \ amg_d_l1_jac_smoother_clone.o \ +amg_d_poly_smoother_apply_vect.o \ +amg_d_poly_smoother_bld.o \ +amg_d_poly_smoother_cnv.o \ +amg_d_poly_smoother_clone.o \ +amg_d_poly_smoother_clone_settings.o \ +amg_d_poly_smoother_clear_data.o \ +amg_d_poly_smoother_descr.o \ +amg_d_poly_smoother_dmp.o \ +amg_d_poly_smoother_csetc.o \ +amg_d_poly_smoother_cseti.o \ +amg_d_poly_smoother_csetr.o \ amg_s_as_smoother_apply.o \ amg_s_as_smoother_apply_vect.o \ amg_s_as_smoother_bld.o \ diff --git a/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 b/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 index 21b6da06..1a084f87 100644 --- a/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 +++ b/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 @@ -72,6 +72,7 @@ subroutine amg_d_jac_smoother_clone(sm,smout,info) smo%checkiter = sm%checkiter smo%printiter = sm%printiter smo%tol = sm%tol + smo%pa => sm%pa call sm%nd%clone(smo%nd,info) if ((info==psb_success_).and.(allocated(sm%sv))) then allocate(smout%sv,mold=sm%sv,stat=info) diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 new file mode 100644 index 00000000..2fc3115c --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 @@ -0,0 +1,357 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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. +! +! +subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& + & sweeps,work,wv,info,init,initu) + + use psb_base_mod + use amg_d_diag_solver + use psb_base_krylov_conv_mod, only : log_conv + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_d_poly_smoother_type), intent(inout) :: sm + 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 + integer(psb_ipk_), intent(in) :: sweeps + 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 + ! + integer(psb_ipk_) :: n_row,n_col + type(psb_d_vect_type) :: tx, ty, tz, r + real(psb_dpk_), pointer :: aux(:) + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: np, me, i, err_act + character :: trans_, init_ + real(psb_dpk_) :: res, resdenum + real(psb_dpk_) :: cz, cr + character(len=20) :: name='d_poly_smoother_apply_v' + + call psb_erractionsave(err_act) + + info = psb_success_ + ctxt = desc_data%get_context() + call psb_info(ctxt,me,np) + + + if (present(init)) then + init_ = psb_toupper(init) + else + init_='Z' + end if + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (4*n_col <= size(work)) then + aux => work(:) + else + allocate(aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/4*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + endif +!!$ if (me == 0) write(0,*) name,' Unimplemented apply_vect ' +!!$ info =psb_err_internal_error_ +!!$ call psb_errpush(info,& +!!$ & name,a_err='Error in sub_aply Polynomial not implemented') +!!$ goto 9999 + + if (size(wv) < 4) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + + associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4)) + + call psb_geaxpby(done,x,dzero,r,desc_data,info) + call tx%zero() + call ty%zero() + call tz%zero() + + do i=1, sm%pdegree + call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Y') + cz = (2*i*done-3)/(2*i*done+done) + cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) + call psb_geaxpby(cr,ty,cz,tz,desc_data,info) + call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info) + call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_) + end do + + if (info == psb_success_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info) + + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='polynomial smoother') + goto 9999 + end if + end associate + + + + +!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then +!!$ ! if .not.sv%is_iterative, there's no need to pass init +!!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') +!!$ goto 9999 +!!$ endif +!!$ +!!$ else if (sweeps >= 0) then +!!$ select type (smsv => sm%sv) +!!$ class is (amg_d_diag_solver_type) +!!$ ! +!!$ ! This means we are dealing with a pure Jacobi smoother/solver. +!!$ ! +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,y,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,initu,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), +!!$ ! where is the diagonal and A the matrix. +!!$ ! +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(done,x,dzero,r,r,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if ( res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ +!!$ end associate +!!$ +!!$ class default +!!$ ! +!!$ ! +!!$ ! Apply multiple sweeps of a block-Jacobi solver +!!$ ! to compute an approximate solution of a linear system. +!!$ ! +!!$ ! +!!$ if (size(wv) < 2) then +!!$ info = psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='invalid wv size in smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ +!!$ ! +!!$ ! Unroll the first iteration and fold it inside SELECT CASE +!!$ ! this will save one AXPBY and one SPMM when INIT=Z, and will be +!!$ ! significant when sweeps=1 (a common case) +!!$ ! +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,y,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_geaxpby(done,initu,dzero,ty,desc_data,info) +!!$ call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the +!!$ ! block diagonal part and the remaining part of the local matrix +!!$ ! and Y(j) is the approximate solution at sweep j. +!!$ ! +!!$ call psb_geaxpby(done,x,dzero,tx,desc_data,info) +!!$ call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(done,x,dzero,r,r,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if (res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ end associate +!!$ end select +!!$ +!!$ else +!!$ +!!$ info = psb_err_iarg_neg_ +!!$ call psb_errpush(info,name,& +!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) +!!$ goto 9999 +!!$ +!!$ endif +!!$ + if (.not.(4*n_col <= size(work))) then + deallocate(aux) + endif + +!!$ if(sm%checkres) then +!!$ call psb_gefree(r,desc_data,info) +!!$ end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + + end subroutine amg_d_poly_smoother_apply_vect diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 new file mode 100644 index 00000000..be83a0ae --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 @@ -0,0 +1,107 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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. +! +! +subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) + + use psb_base_mod + use amg_d_diag_solver + use amg_d_beta_coeff_mod + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_bld + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_d_poly_smoother_type), intent(inout) :: sm + 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 + ! Local variables + type(psb_dspmat_type) :: tmpa + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level + character(len=20) :: name='d_poly_smoother_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' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + call psb_realloc(sm%pdegree,sm%poly_beta,info) + sm%poly_beta(1:sm%pdegree) = amg_d_beta_mat(1:sm%pdegree,sm%pdegree) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%degree') + goto 9999 + end if + sm%pa => a + call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='sv%build') + goto 9999 + end if + + if (sm%rho_ba <= dzero) then + sm%rho_ba = psb_dspnrm1(a,desc_a,info) + 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_poly_smoother_bld diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 new file mode 100644 index 00000000..ac526bca --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 @@ -0,0 +1,70 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! 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. +! +! +subroutine amg_d_poly_smoother_clear_data(sm,info) + + use psb_base_mod + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_clear_data + Implicit None + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_poly_smoother_clear_data' + + call psb_erractionsave(err_act) + + info = 0 + sm%pdegree = 0 + if (allocated(sm%poly_beta)) deallocate(sm%poly_beta) + sm%pa => null() + if ((info==0).and.allocated(sm%sv)) then + call sm%sv%clear_data(info) + end if + if (info /= 0) then + info = psb_err_internal_error_ + 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 amg_d_poly_smoother_clear_data diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_clone.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_clone.f90 new file mode 100644 index 00000000..7dfb2655 --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_clone.f90 @@ -0,0 +1,90 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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. +! +! +subroutine amg_d_poly_smoother_clone(sm,smout,info) + + use psb_base_mod + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_clone + + Implicit None + + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + class(amg_d_base_smoother_type), allocatable, intent(inout) :: smout + integer(psb_ipk_), intent(out) :: info + ! Local variables + integer(psb_ipk_) :: err_act + + + info=psb_success_ + call psb_erractionsave(err_act) + + if (allocated(smout)) then + call smout%free(info) + if (info == psb_success_) deallocate(smout, stat=info) + end if + if (info == psb_success_) & + & allocate(amg_d_poly_smoother_type :: smout, stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + goto 9999 + end if + + select type(smo => smout) + type is (amg_d_poly_smoother_type) + smo%pdegree = sm%pdegree + smo%rho_ba = sm%rho_ba + smo%poly_beta = sm%poly_beta + smo%pa => sm%pa + if ((info==psb_success_).and.(allocated(sm%sv))) then + allocate(smout%sv,mold=sm%sv,stat=info) + if (info == psb_success_) call sm%sv%clone(smo%sv,info) + end if + + class default + info = psb_err_internal_error_ + end select + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_d_poly_smoother_clone diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 new file mode 100644 index 00000000..1fbdac37 --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 @@ -0,0 +1,96 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! asd on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! 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. +! +! +subroutine amg_d_poly_smoother_clone_settings(sm,smout,info) + + use psb_base_mod + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_clone_settings + Implicit None + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + class(amg_d_base_smoother_type), intent(inout) :: smout + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_poly_smoother_clone_settings' + + call psb_erractionsave(err_act) + + info = psb_success_ + + select type(smout) + class is(amg_d_poly_smoother_type) + + smout%pa => null() + smout%pdegree = sm%pdegree + + if (allocated(smout%sv)) then + if (.not.same_type_as(sm%sv,smout%sv)) then + call smout%sv%free(info) + if (info == 0) deallocate(smout%sv,stat=info) + end if + end if + if (info /= 0) then + info = psb_err_internal_error_ + else + if (allocated(smout%sv)) then + if (same_type_as(sm%sv,smout%sv)) then + call sm%sv%clone_settings(smout%sv,info) + else + info = psb_err_internal_error_ + end if + else + allocate(smout%sv,mold=sm%sv,stat=info) + if (info == 0) call sm%sv%clone_settings(smout%sv,info) + if (info /= 0) info = psb_err_internal_error_ + end if + end if + class default + info = psb_err_internal_error_ + end select + + if (info /= 0) then + 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 amg_d_poly_smoother_clone_settings diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_cnv.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_cnv.f90 new file mode 100644 index 00000000..42d2ce15 --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_cnv.f90 @@ -0,0 +1,77 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! 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. +! +! +subroutine amg_d_poly_smoother_cnv(sm,info,amold,vmold,imold) + + use psb_base_mod + use amg_d_diag_solver + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_cnv + Implicit None + + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + 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 + ! Local variables + integer(psb_ipk_) :: i, err_act, debug_unit, debug_level + character(len=20) :: name='d_poly_smoother_cnv', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + + if (allocated(sm%sv)) & + & call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='solver cnv') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine amg_d_poly_smoother_cnv diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 new file mode 100644 index 00000000..3c1fef00 --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 @@ -0,0 +1,72 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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. +! +! +subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx) + + use psb_base_mod + use amg_d_poly_smoother, amg_protect_nam => amg_d_poly_smoother_csetc + Implicit None + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + 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_poly_smoother_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(trim(what))) + case default + call sm%amg_d_base_smoother_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 amg_d_poly_smoother_csetc diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 new file mode 100644 index 00000000..8bc48724 --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 @@ -0,0 +1,69 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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. +! +! +subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx) + + use psb_base_mod + use amg_d_poly_smoother, amg_protect_nam => amg_d_poly_smoother_cseti + Implicit None + + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + 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_poly_smoother_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + case('SMOOTHER_DEGREE') + sm%pdegree = val + case default + call sm%amg_d_base_smoother_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 amg_d_poly_smoother_cseti diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 new file mode 100644 index 00000000..de308a8e --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90 @@ -0,0 +1,69 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! 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. +! +! +subroutine amg_d_poly_smoother_csetr(sm,what,val,info,idx) + + use psb_base_mod + use amg_d_poly_smoother, amg_protect_nam => amg_d_poly_smoother_csetr + Implicit None + + ! Arguments + class(amg_d_poly_smoother_type), intent(inout) :: sm + 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_poly_smoother_csetr' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + case('RHO_BA') + sm%rho_ba = val + case default + call sm%amg_d_base_smoother_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 amg_d_poly_smoother_csetr diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 new file mode 100644 index 00000000..d17c5aa5 --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 @@ -0,0 +1,95 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! 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. +! +! +subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) + + use psb_base_mod + use amg_d_diag_solver + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_descr + use amg_d_diag_solver + use amg_d_gs_solver + + Implicit None + + ! Arguments + class(amg_d_poly_smoother_type), intent(in) :: sm + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + character(len=*), intent(in), optional :: prefix + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20), parameter :: name='amg_d_poly_smoother_descr' + integer(psb_ipk_) :: iout_ + logical :: coarse_ + character(1024) :: prefix_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(coarse)) then + coarse_ = coarse + else + coarse_ = .false. + end if + if (present(iout)) then + iout_ = iout + else + iout_ = psb_out_unit + endif + if (present(prefix)) then + prefix_ = prefix + else + prefix_ = "" + end if + + write(iout_,*) trim(prefix_), ' Polynomial smoother ' + write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree + write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba + if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), ' Coefficients: ',sm%poly_beta(1:sm%pdegree) + if (allocated(sm%sv)) then + write(iout_,*) trim(prefix_), ' Local solver details:' + call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix) + end if + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_d_poly_smoother_descr diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_dmp.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_dmp.f90 new file mode 100644 index 00000000..19144f07 --- /dev/null +++ b/amgprec/impl/smoother/amg_d_poly_smoother_dmp.f90 @@ -0,0 +1,90 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! 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. +! +! +subroutine amg_d_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solver,global_num) + + use psb_base_mod + use amg_d_poly_smoother, amg_protect_nam => amg_d_poly_smoother_dmp + implicit none + class(amg_d_poly_smoother_type), intent(in) :: sm + 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) :: smoother, solver, global_num + integer(psb_ipk_) :: i, j, il1, iln, lname, lev + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: iam, np + character(len=80) :: prefix_ + character(len=120) :: fname ! len should be at least 20 more than + integer(psb_lpk_), allocatable :: iv(:) + logical :: smoother_, global_num_ + ! len of prefix_ + + info = 0 + + if (present(prefix)) then + prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) + else + prefix_ = "dump_smth_d" + end if + ctxt = desc%get_context() + call psb_info(ctxt,iam,np) + + if (present(smoother)) then + smoother_ = smoother + else + smoother_ = .false. + end if + if (present(global_num)) then + global_num_ = global_num + else + global_num_ = .false. + end if + lname = len_trim(prefix_) + fname = trim(prefix_) + write(fname(lname+1:lname+5),'(a,i3.3)') '_poly',iam + lname = lname + 8 + ! to be completed + + + + ! At base level do nothing for the smoother + if (allocated(sm%sv)) & + & call sm%sv%dump(desc,level,info,solver=solver,prefix=prefix,global_num=global_num) + +end subroutine amg_d_poly_smoother_dmp diff --git a/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 index 388afe1e..9f6109db 100644 --- a/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 @@ -52,7 +52,7 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) class(psb_c_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables - integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype !!$ complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level @@ -97,10 +97,24 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 endif - + select case(sv%fact_type) + case (amg_ilu_n_) + psb_fctype = amg_ilu_n_ + case (amg_milu_n_) + psb_fctype = amg_milu_n_ + case (amg_ilu_t_) + psb_fctype = amg_ilu_t_ + case default + ! If we end up here, something was wrong up in the call chain. + info = psb_err_input_value_invalid_i_ + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ithree,sv%fact_type,izero,izero,izero/)) + goto 9999 + end select + select case(sv%fact_type) - case (psb_ilu_t_) + case (amg_ilu_t_) ! ! ILU(k,t) ! @@ -124,7 +138,7 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 end if - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) ! ! ILU(k) and MILU(k) ! @@ -140,17 +154,17 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! There seems to be a problem with the separate implementation of MILU(0), ! contained into psb_ilu0_fact. This must be investigated. For the time being, ! resort to the implementation of MILU(k) with k=0. - if (sv%fact_type == psb_ilu_n_) then - call psb_ilu0_fact(sv%fact_type,a,sv%l,sv%u,& + if (sv%fact_type == amg_ilu_n_) then + call psb_ilu0_fact(psb_fctype,a,sv%l,sv%u,& & sv%d,info,blck=b) else - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) endif case(1:) ! Fill-in >= 1 ! The same routine implements both ILU(k) and MILU(k) - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) end select if (info /= psb_success_) then diff --git a/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 index 510ea26a..c4ab246a 100644 --- a/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 @@ -52,7 +52,7 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) class(psb_d_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables - integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype !!$ real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level @@ -97,10 +97,24 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 endif - + select case(sv%fact_type) + case (amg_ilu_n_) + psb_fctype = amg_ilu_n_ + case (amg_milu_n_) + psb_fctype = amg_milu_n_ + case (amg_ilu_t_) + psb_fctype = amg_ilu_t_ + case default + ! If we end up here, something was wrong up in the call chain. + info = psb_err_input_value_invalid_i_ + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ithree,sv%fact_type,izero,izero,izero/)) + goto 9999 + end select + select case(sv%fact_type) - case (psb_ilu_t_) + case (amg_ilu_t_) ! ! ILU(k,t) ! @@ -124,7 +138,7 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 end if - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) ! ! ILU(k) and MILU(k) ! @@ -140,17 +154,17 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! There seems to be a problem with the separate implementation of MILU(0), ! contained into psb_ilu0_fact. This must be investigated. For the time being, ! resort to the implementation of MILU(k) with k=0. - if (sv%fact_type == psb_ilu_n_) then - call psb_ilu0_fact(sv%fact_type,a,sv%l,sv%u,& + if (sv%fact_type == amg_ilu_n_) then + call psb_ilu0_fact(psb_fctype,a,sv%l,sv%u,& & sv%d,info,blck=b) else - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) endif case(1:) ! Fill-in >= 1 ! The same routine implements both ILU(k) and MILU(k) - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) end select if (info /= psb_success_) then diff --git a/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 index 4516b134..8cfadb3b 100644 --- a/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 @@ -52,7 +52,7 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) class(psb_s_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables - integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype !!$ real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level @@ -97,10 +97,24 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 endif - + select case(sv%fact_type) + case (amg_ilu_n_) + psb_fctype = amg_ilu_n_ + case (amg_milu_n_) + psb_fctype = amg_milu_n_ + case (amg_ilu_t_) + psb_fctype = amg_ilu_t_ + case default + ! If we end up here, something was wrong up in the call chain. + info = psb_err_input_value_invalid_i_ + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ithree,sv%fact_type,izero,izero,izero/)) + goto 9999 + end select + select case(sv%fact_type) - case (psb_ilu_t_) + case (amg_ilu_t_) ! ! ILU(k,t) ! @@ -124,7 +138,7 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 end if - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) ! ! ILU(k) and MILU(k) ! @@ -140,17 +154,17 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! There seems to be a problem with the separate implementation of MILU(0), ! contained into psb_ilu0_fact. This must be investigated. For the time being, ! resort to the implementation of MILU(k) with k=0. - if (sv%fact_type == psb_ilu_n_) then - call psb_ilu0_fact(sv%fact_type,a,sv%l,sv%u,& + if (sv%fact_type == amg_ilu_n_) then + call psb_ilu0_fact(psb_fctype,a,sv%l,sv%u,& & sv%d,info,blck=b) else - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) endif case(1:) ! Fill-in >= 1 ! The same routine implements both ILU(k) and MILU(k) - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) end select if (info /= psb_success_) then diff --git a/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 index aead8d16..3afa82f7 100644 --- a/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 @@ -52,7 +52,7 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) class(psb_z_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables - integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype !!$ complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level @@ -97,10 +97,24 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 endif - + select case(sv%fact_type) + case (amg_ilu_n_) + psb_fctype = amg_ilu_n_ + case (amg_milu_n_) + psb_fctype = amg_milu_n_ + case (amg_ilu_t_) + psb_fctype = amg_ilu_t_ + case default + ! If we end up here, something was wrong up in the call chain. + info = psb_err_input_value_invalid_i_ + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ithree,sv%fact_type,izero,izero,izero/)) + goto 9999 + end select + select case(sv%fact_type) - case (psb_ilu_t_) + case (amg_ilu_t_) ! ! ILU(k,t) ! @@ -124,7 +138,7 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) goto 9999 end if - case(psb_ilu_n_,psb_milu_n_) + case(amg_ilu_n_,amg_milu_n_) ! ! ILU(k) and MILU(k) ! @@ -140,17 +154,17 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! There seems to be a problem with the separate implementation of MILU(0), ! contained into psb_ilu0_fact. This must be investigated. For the time being, ! resort to the implementation of MILU(k) with k=0. - if (sv%fact_type == psb_ilu_n_) then - call psb_ilu0_fact(sv%fact_type,a,sv%l,sv%u,& + if (sv%fact_type == amg_ilu_n_) then + call psb_ilu0_fact(psb_fctype,a,sv%l,sv%u,& & sv%d,info,blck=b) else - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) endif case(1:) ! Fill-in >= 1 ! The same routine implements both ILU(k) and MILU(k) - call psb_iluk_fact(sv%fill_in,sv%fact_type,& + call psb_iluk_fact(sv%fill_in,psb_fctype,& & a,sv%l,sv%u,sv%d,info,blck=b) end select if (info /= psb_success_) then diff --git a/samples/advanced/pdegen/amg_d_pde3d.F90 b/samples/advanced/pdegen/amg_d_pde3d.F90 index 006c6d6f..410b6e01 100644 --- a/samples/advanced/pdegen/amg_d_pde3d.F90 +++ b/samples/advanced/pdegen/amg_d_pde3d.F90 @@ -144,6 +144,7 @@ program amg_d_pde3d ! AMG smoother or pre-smoother; also 1-lev preconditioner character(len=16) :: smther ! (pre-)smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps + integer(psb_ipk_) :: degree ! degree for polynomial smoother integer(psb_ipk_) :: novr ! number of overlap layers character(len=16) :: restr ! restriction over application of AS character(len=16) :: prol ! prolongation over application of AS @@ -158,6 +159,7 @@ program amg_d_pde3d ! AMG post-smoother; ignored by 1-lev preconditioner character(len=16) :: smther2 ! post-smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps + integer(psb_ipk_) :: degree2 ! degree for polynomial smoother integer(psb_ipk_) :: novr2 ! number of overlap layers character(len=16) :: restr2 ! restriction over application of AS character(len=16) :: prol2 ! prolongation over application of AS @@ -285,10 +287,11 @@ program amg_d_pde3d ! 1-level sweeps from "outer_sweeps" call prec%set('smoother_sweeps', p_choice%jsweeps, info) - case ('BJAC') + case ('BJAC','POLY') call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('sub_solve', p_choice%solve, info) call prec%set('solver_sweeps', p_choice%ssweeps, info) + call prec%set('smoother_degree', p_choice%degree, info) if (psb_toupper(p_choice%solve)=='MUMPS') & & call prec%set('mumps_loc_glob','local_solver',info) call prec%set('sub_fillin', p_choice%fill, info) @@ -336,7 +339,8 @@ program amg_d_pde3d call prec%set('smoother_type', p_choice%smther, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + call prec%set('smoother_degree', p_choice%degree, info) + select case (psb_toupper(p_choice%smther)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -366,6 +370,7 @@ program amg_d_pde3d if (psb_toupper(p_choice%smther2) /= 'NONE') then call prec%set('smoother_type', p_choice%smther2, info,pos='post') call prec%set('smoother_sweeps', p_choice%jsweeps2, info,pos='post') + call prec%set('smoother_degree', p_choice%degree2, info,pos='post') select case (psb_toupper(p_choice%smther2)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -402,6 +407,7 @@ program amg_d_pde3d call prec%set('coarse_sweeps', p_choice%cjswp, info) end select +!!$ call prec%descr(info,iout=psb_out_unit) ! build the preconditioner call psb_barrier(ctxt) @@ -581,6 +587,7 @@ contains ! First smoother / 1-lev preconditioner call read_data(prec%smther,inp_unit) ! smoother type call read_data(prec%jsweeps,inp_unit) ! (pre-)smoother / 1-lev prec sweeps + call read_data(prec%degree,inp_unit) ! (pre-)smoother / 1-lev prec sweeps call read_data(prec%novr,inp_unit) ! number of overlap layers call read_data(prec%restr,inp_unit) ! restriction over application of AS call read_data(prec%prol,inp_unit) ! prolongation over application of AS @@ -593,11 +600,12 @@ contains ! Second smoother/ AMG post-smoother (if NONE ignored in main) call read_data(prec%smther2,inp_unit) ! smoother type call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps + call read_data(prec%degree2,inp_unit) ! (post-)smoother sweeps call read_data(prec%novr2,inp_unit) ! number of overlap layers call read_data(prec%restr2,inp_unit) ! restriction over application of AS call read_data(prec%prol2,inp_unit) ! prolongation over application of AS call read_data(prec%solve2,inp_unit) ! local subsolver - call read_data(prec%ssweeps2,inp_unit) ! inner solver sweeps + call read_data(prec%ssweeps2,inp_unit) ! inner solver sweeps call read_data(prec%variant2,inp_unit) ! AINV variant call read_data(prec%fill2,inp_unit) ! fill-in for incomplete LU call read_data(prec%invfill2,inp_unit) !Inverse fill-in for INVK @@ -663,6 +671,7 @@ contains ! broadcast first (pre-)smoother / 1-lev prec data call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) + call psb_bcast(ctxt,prec%degree) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) call psb_bcast(ctxt,prec%prol) @@ -675,6 +684,7 @@ contains ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) + call psb_bcast(ctxt,prec%degree2) call psb_bcast(ctxt,prec%novr2) call psb_bcast(ctxt,prec%restr2) call psb_bcast(ctxt,prec%prol2)