From b331b1b9280ae8f6379c0553caf61b743423191c Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 4 May 2018 14:36:23 +0100 Subject: [PATCH] Merge aggregator: module interfaces now compile. --- mlprec/Makefile | 36 +++- mlprec/mld_base_prec_type.F90 | 13 +- mlprec/mld_c_base_aggregator_mod.f90 | 209 +++++++++++++++++++ mlprec/mld_c_hybrid_aggregator_mod.F90 | 131 ++++++++++++ mlprec/mld_c_inner_mod.f90 | 25 +++ mlprec/mld_c_onelev_mod.f90 | 49 ++++- mlprec/mld_c_prec_mod.f90 | 11 +- mlprec/mld_c_prec_type.f90 | 17 +- mlprec/mld_c_symdec_aggregator_mod.f90 | 128 ++++++++++++ mlprec/mld_d_base_aggregator_mod.f90 | 209 +++++++++++++++++++ mlprec/mld_d_bcmatch_aggregator_mod.F90 | 262 ++++++++++++++++++++++++ mlprec/mld_d_hybrid_aggregator_mod.F90 | 131 ++++++++++++ mlprec/mld_d_inner_mod.f90 | 25 +++ mlprec/mld_d_onelev_mod.f90 | 49 ++++- mlprec/mld_d_prec_mod.f90 | 11 +- mlprec/mld_d_prec_type.f90 | 17 +- mlprec/mld_d_symdec_aggregator_mod.f90 | 128 ++++++++++++ mlprec/mld_s_base_aggregator_mod.f90 | 209 +++++++++++++++++++ mlprec/mld_s_hybrid_aggregator_mod.F90 | 131 ++++++++++++ mlprec/mld_s_inner_mod.f90 | 25 +++ mlprec/mld_s_onelev_mod.f90 | 49 ++++- mlprec/mld_s_prec_mod.f90 | 11 +- mlprec/mld_s_prec_type.f90 | 17 +- mlprec/mld_s_symdec_aggregator_mod.f90 | 128 ++++++++++++ mlprec/mld_z_base_aggregator_mod.f90 | 209 +++++++++++++++++++ mlprec/mld_z_hybrid_aggregator_mod.F90 | 131 ++++++++++++ mlprec/mld_z_inner_mod.f90 | 25 +++ mlprec/mld_z_onelev_mod.f90 | 49 ++++- mlprec/mld_z_prec_mod.f90 | 11 +- mlprec/mld_z_prec_type.f90 | 17 +- mlprec/mld_z_symdec_aggregator_mod.f90 | 128 ++++++++++++ 31 files changed, 2523 insertions(+), 68 deletions(-) create mode 100644 mlprec/mld_c_base_aggregator_mod.f90 create mode 100644 mlprec/mld_c_hybrid_aggregator_mod.F90 create mode 100644 mlprec/mld_c_symdec_aggregator_mod.f90 create mode 100644 mlprec/mld_d_base_aggregator_mod.f90 create mode 100644 mlprec/mld_d_bcmatch_aggregator_mod.F90 create mode 100644 mlprec/mld_d_hybrid_aggregator_mod.F90 create mode 100644 mlprec/mld_d_symdec_aggregator_mod.f90 create mode 100644 mlprec/mld_s_base_aggregator_mod.f90 create mode 100644 mlprec/mld_s_hybrid_aggregator_mod.F90 create mode 100644 mlprec/mld_s_symdec_aggregator_mod.f90 create mode 100644 mlprec/mld_z_base_aggregator_mod.f90 create mode 100644 mlprec/mld_z_hybrid_aggregator_mod.F90 create mode 100644 mlprec/mld_z_symdec_aggregator_mod.f90 diff --git a/mlprec/Makefile b/mlprec/Makefile index e853352a..d89347b7 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -11,25 +11,33 @@ DMODOBJS=mld_d_prec_type.o mld_d_ilu_fact_mod.o \ mld_d_inner_mod.o mld_d_ilu_solver.o mld_d_diag_solver.o mld_d_jac_smoother.o mld_d_as_smoother.o \ mld_d_umf_solver.o mld_d_slu_solver.o mld_d_sludist_solver.o mld_d_id_solver.o\ mld_d_base_solver_mod.o mld_d_base_smoother_mod.o mld_d_onelev_mod.o \ - mld_d_gs_solver.o mld_d_mumps_solver.o + mld_d_gs_solver.o mld_d_mumps_solver.o \ + mld_d_base_aggregator_mod.o mld_d_hybrid_aggregator_mod.o \ + mld_d_symdec_aggregator_mod.o mld_d_bcmatch_aggregator_mod.o SMODOBJS=mld_s_prec_type.o mld_s_ilu_fact_mod.o \ mld_s_inner_mod.o mld_s_ilu_solver.o mld_s_diag_solver.o mld_s_jac_smoother.o mld_s_as_smoother.o \ mld_s_slu_solver.o mld_s_id_solver.o\ mld_s_base_solver_mod.o mld_s_base_smoother_mod.o mld_s_onelev_mod.o \ - mld_s_gs_solver.o mld_s_mumps_solver.o + mld_s_gs_solver.o mld_s_mumps_solver.o \ + mld_s_base_aggregator_mod.o mld_s_hybrid_aggregator_mod.o \ + mld_s_symdec_aggregator_mod.o ZMODOBJS=mld_z_prec_type.o mld_z_ilu_fact_mod.o \ mld_z_inner_mod.o mld_z_ilu_solver.o mld_z_diag_solver.o mld_z_jac_smoother.o mld_z_as_smoother.o \ mld_z_umf_solver.o mld_z_slu_solver.o mld_z_sludist_solver.o mld_z_id_solver.o\ mld_z_base_solver_mod.o mld_z_base_smoother_mod.o mld_z_onelev_mod.o \ - mld_z_gs_solver.o mld_z_mumps_solver.o + mld_z_gs_solver.o mld_z_mumps_solver.o \ + mld_z_base_aggregator_mod.o mld_z_hybrid_aggregator_mod.o \ + mld_z_symdec_aggregator_mod.o CMODOBJS=mld_c_prec_type.o mld_c_ilu_fact_mod.o \ mld_c_inner_mod.o mld_c_ilu_solver.o mld_c_diag_solver.o mld_c_jac_smoother.o mld_c_as_smoother.o \ mld_c_slu_solver.o mld_c_id_solver.o\ mld_c_base_solver_mod.o mld_c_base_smoother_mod.o mld_c_onelev_mod.o \ - mld_c_gs_solver.o mld_c_mumps_solver.o + mld_c_gs_solver.o mld_c_mumps_solver.o \ + mld_c_base_aggregator_mod.o mld_c_hybrid_aggregator_mod.o \ + mld_c_symdec_aggregator_mod.o @@ -84,10 +92,22 @@ mld_d_prec_type.o: mld_d_onelev_mod.o mld_c_prec_type.o: mld_c_onelev_mod.o mld_z_prec_type.o: mld_z_onelev_mod.o -mld_s_onelev_mod.o: mld_s_base_smoother_mod.o -mld_d_onelev_mod.o: mld_d_base_smoother_mod.o -mld_c_onelev_mod.o: mld_c_base_smoother_mod.o -mld_z_onelev_mod.o: mld_z_base_smoother_mod.o +mld_s_onelev_mod.o: mld_s_base_smoother_mod.o mld_s_base_aggregator_mod.o +mld_d_onelev_mod.o: mld_d_base_smoother_mod.o mld_d_base_aggregator_mod.o +mld_c_onelev_mod.o: mld_c_base_smoother_mod.o mld_c_base_aggregator_mod.o +mld_z_onelev_mod.o: mld_z_base_smoother_mod.o mld_z_base_aggregator_mod.o + +mld_s_base_aggregator_mod.o: mld_base_prec_type.o +mld_s_hybrid_aggregator_mod.o mld_s_symdec_aggregator_mod.o: mld_s_base_aggregator_mod.o + +mld_d_base_aggregator_mod.o: mld_base_prec_type.o +mld_d_bcmatch_aggregator_mod.o mld_d_hybrid_aggregator_mod.o mld_d_symdec_aggregator_mod.o: mld_d_base_aggregator_mod.o + +mld_c_base_aggregator_mod.o: mld_base_prec_type.o +mld_c_hybrid_aggregator_mod.o mld_c_symdec_aggregator_mod.o: mld_c_base_aggregator_mod.o + +mld_z_base_aggregator_mod.o: mld_base_prec_type.o +mld_z_hybrid_aggregator_mod.o mld_z_symdec_aggregator_mod.o: mld_z_base_aggregator_mod.o mld_s_base_smoother_mod.o: mld_s_base_solver_mod.o mld_d_base_smoother_mod.o: mld_d_base_solver_mod.o diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index 5de5c78c..d473cb2d 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -534,7 +534,7 @@ contains ! Routines printing out a description of the preconditioner ! - subroutine ml_parms_mldescr(pm,iout,info) + subroutine ml_parms_mldescr(pm,iout,info,aggr_name) Implicit None @@ -542,7 +542,7 @@ contains class(mld_ml_parms), intent(in) :: pm integer(psb_ipk_), intent(in) :: iout integer(psb_ipk_), intent(out) :: info - + character(len=*), intent(in), optional :: aggr_name info = psb_success_ if ((pm%ml_cycle>=mld_no_ml_).and.(pm%ml_cycle<=mld_max_ml_cycle_)) then @@ -558,8 +558,13 @@ contains & pm%sweeps_pre ,' post: ', pm%sweeps_post end select - write(iout,*) ' Aggregation type: ',& - & aggr_type_names(pm%aggr_type) + if (present(aggr_name)) then + write(iout,*) ' Aggregation type: ', & + & aggr_name + else + write(iout,*) ' Aggregation type: ',& + & aggr_type_names(pm%aggr_type) + end if write(iout,*) ' parallel algorithm: ',& & par_aggr_alg_names(pm%par_aggr_alg) if (pm%par_aggr_alg /= mld_ext_aggr_) then diff --git a/mlprec/mld_c_base_aggregator_mod.f90 b/mlprec/mld_c_base_aggregator_mod.f90 new file mode 100644 index 00000000..6fd291cc --- /dev/null +++ b/mlprec/mld_c_base_aggregator_mod.f90 @@ -0,0 +1,209 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. The basic version is the +! decoupled aggregation algorithm presented in +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_c_base_aggregator_mod + + use mld_base_prec_type, only : mld_sml_parms + use psb_base_mod, only : psb_cspmat_type, psb_c_vect_type, & + & psb_c_base_vect_type, psb_clinmap_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & + & psb_erractionsave, psb_error_handler, psb_success_ + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type mld_c_base_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_c_base_aggregator_build_tprol + procedure, pass(ag) :: mat_asb => mld_c_base_aggregator_mat_asb + procedure, pass(ag) :: update_level => mld_c_base_aggregator_update_level + procedure, pass(ag) :: clone => mld_c_base_aggregator_clone + procedure, pass(ag) :: free => mld_c_base_aggregator_free + procedure, pass(ag) :: default => mld_c_base_aggregator_default + procedure, pass(ag) :: descr => mld_c_base_aggregator_descr + procedure, nopass :: fmt => mld_c_base_aggregator_fmt + end type mld_c_base_aggregator_type + + + interface + subroutine mld_c_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_c_base_aggregator_type, psb_desc_type, psb_cspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_c_base_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_base_aggregator_build_tprol + end interface + + interface + subroutine mld_c_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,& + & op_prol,op_restr,info) + import :: mld_c_base_aggregator_type, psb_desc_type, psb_cspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_c_base_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(psb_cspmat_type), intent(inout) :: op_prol + type(psb_cspmat_type), intent(out) :: ac,op_restr + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_base_aggregator_mat_asb + end interface + +contains + + subroutine mld_c_base_aggregator_update_level(ag,agnext,info) + implicit none + class(mld_c_base_aggregator_type), target, intent(inout) :: ag, agnext + integer(psb_ipk_), intent(out) :: info + + ! + ! Base version does nothing. + ! + info = 0 + end subroutine mld_c_base_aggregator_update_level + + subroutine mld_c_base_aggregator_clone(ag,agnext,info) + implicit none + class(mld_c_base_aggregator_type), intent(inout) :: ag + class(mld_c_base_aggregator_type), allocatable, intent(inout) :: agnext + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(agnext)) then + call agnext%free(info) + if (info == 0) deallocate(agnext,stat=info) + end if + if (info /= 0) return + allocate(agnext,source=ag,stat=info) + + end subroutine mld_c_base_aggregator_clone + + subroutine mld_c_base_aggregator_free(ag,info) + implicit none + class(mld_c_base_aggregator_type), intent(inout) :: ag + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + return + end subroutine mld_c_base_aggregator_free + + subroutine mld_c_base_aggregator_default(ag) + implicit none + class(mld_c_base_aggregator_type), intent(inout) :: ag + + ! Here we need do nothing + + return + end subroutine mld_c_base_aggregator_default + + function mld_c_base_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Decoupled aggregation" + end function mld_c_base_aggregator_fmt + + subroutine mld_c_base_aggregator_descr(ag,parms,iout,info) + implicit none + class(mld_c_base_aggregator_type), intent(in) :: ag + type(mld_sml_parms), intent(in) :: parms + integer(psb_ipk_), intent(in) :: iout + integer(psb_ipk_), intent(out) :: info + + call parms%mldescr(iout,info,aggr_name=ag%fmt()) + + return + end subroutine mld_c_base_aggregator_descr + +end module mld_c_base_aggregator_mod diff --git a/mlprec/mld_c_hybrid_aggregator_mod.F90 b/mlprec/mld_c_hybrid_aggregator_mod.F90 new file mode 100644 index 00000000..bbc62953 --- /dev/null +++ b/mlprec/mld_c_hybrid_aggregator_mod.F90 @@ -0,0 +1,131 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This variant is based on the hybrid method +! presented in +! +! S. Gratton, P. Henon, P. Jiranek and X. Vasseur: +! Reducing complexity of algebraic multigrid by aggregation +! Numerical Lin. Algebra with Applications, 2016, 23:501-518 +! +module mld_c_hybrid_aggregator_mod + + use mld_c_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_c_base_aggregator_type) :: mld_c_hybrid_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_c_hybrid_aggregator_build_tprol +!!$ procedure, pass(ag) :: mat_asb => mld_c_base_aggregator_mat_asb +!!$ procedure, pass(ag) :: update_level => mld_c_base_aggregator_update_level +!!$ procedure, pass(ag) :: clone => mld_c_base_aggregator_clone +!!$ procedure, pass(ag) :: free => mld_c_base_aggregator_free +!!$ procedure, pass(ag) :: default => mld_c_base_aggregator_default + procedure, nopass :: fmt => mld_c_hybrid_aggregator_fmt + end type mld_c_hybrid_aggregator_type + + + interface + subroutine mld_c_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_c_hybrid_aggregator_type, psb_desc_type, psb_cspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_c_hybrid_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_hybrid_aggregator_build_tprol + end interface + +contains + + + function mld_c_hybrid_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Hybrid Decoupled aggregation" + end function mld_c_hybrid_aggregator_fmt + + +end module mld_c_hybrid_aggregator_mod diff --git a/mlprec/mld_c_inner_mod.f90 b/mlprec/mld_c_inner_mod.f90 index e5394d4b..0de61dca 100644 --- a/mlprec/mld_c_inner_mod.f90 +++ b/mlprec/mld_c_inner_mod.f90 @@ -137,6 +137,31 @@ module mld_c_inner_mod end subroutine mld_c_dec_map_bld end interface mld_dec_map_bld + interface mld_hyb_map_bld + subroutine mld_c_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: theta + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_hyb_map_bld + end interface mld_hyb_map_bld + + interface mld_map_to_tprol + subroutine mld_c_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + use mld_c_prec_type, only : mld_c_onelev_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_map_to_tprol + end interface mld_map_to_tprol + interface mld_lev_mat_asb subroutine mld_c_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) diff --git a/mlprec/mld_c_onelev_mod.f90 b/mlprec/mld_c_onelev_mod.f90 index f337d4a9..fe5a84a3 100644 --- a/mlprec/mld_c_onelev_mod.f90 +++ b/mlprec/mld_c_onelev_mod.f90 @@ -55,6 +55,7 @@ module mld_c_onelev_mod use mld_base_prec_type use mld_c_base_smoother_mod + use mld_c_base_aggregator_mod use psb_base_mod, only : psb_cspmat_type, psb_c_vect_type, & & psb_c_base_vect_type, psb_clinmap_type, psb_spk_, & & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & @@ -136,9 +137,10 @@ module mld_c_onelev_mod & c_wrk_clone, c_wrk_move_alloc, c_wrk_cnv type mld_c_onelev_type - class(mld_c_base_smoother_type), allocatable :: sm, sm2a - class(mld_c_base_smoother_type), pointer :: sm2 => null() - class(mld_cmlprec_wrk_type), allocatable :: wrk + class(mld_c_base_smoother_type), allocatable :: sm, sm2a + class(mld_c_base_smoother_type), pointer :: sm2 => null() + class(mld_cmlprec_wrk_type), allocatable :: wrk + class(mld_c_base_aggregator_type), allocatable :: aggr type(mld_sml_parms) :: parms type(psb_cspmat_type) :: ac integer(psb_ipk_) :: ac_nz_loc, ac_nz_tot @@ -166,8 +168,9 @@ module mld_c_onelev_mod procedure, pass(lv) :: csetc => mld_c_base_onelev_csetc procedure, pass(lv) :: setsm => mld_c_base_onelev_setsm procedure, pass(lv) :: setsv => mld_c_base_onelev_setsv + procedure, pass(lv) :: setag => mld_c_base_onelev_setag generic, public :: set => seti, setr, setc, & - & cseti, csetr, csetc, setsm, setsv + & cseti, csetr, csetc, setsm, setsv, setag procedure, pass(lv) :: sizeof => c_base_onelev_sizeof procedure, pass(lv) :: get_nzeros => c_base_onelev_get_nzeros procedure, pass(lv) :: get_wrksz => c_base_onelev_get_wrksize @@ -299,6 +302,20 @@ module mld_c_onelev_mod end subroutine mld_c_base_onelev_setsv end interface + interface + subroutine mld_c_base_onelev_setag(lv,val,info,pos) + import :: psb_spk_, mld_c_onelev_type, mld_c_base_aggregator_type, & + & psb_ipk_, psb_long_int_k_, psb_desc_type + Implicit None + + ! Arguments + class(mld_c_onelev_type), target, intent(inout) :: lv + class(mld_c_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + end subroutine mld_c_base_onelev_setag + end interface + interface subroutine mld_c_base_onelev_setc(lv,what,val,info,pos) import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & @@ -323,9 +340,9 @@ module mld_c_onelev_mod class(mld_c_onelev_type), intent(inout) :: lv integer(psb_ipk_), intent(in) :: what - real(psb_spk_), intent(in) :: val + real(psb_spk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info - character(len=*), optional, intent(in) :: pos + character(len=*), optional, intent(in) :: pos end subroutine mld_c_base_onelev_setr end interface @@ -450,7 +467,8 @@ contains Implicit None ! Arguments - class(mld_c_onelev_type), target, intent(inout) :: lv + class(mld_c_onelev_type), target, intent(inout) :: lv + integer(psb_ipk_) :: info lv%parms%sweeps_pre = 1 lv%parms%sweeps_post = 1 @@ -473,7 +491,9 @@ contains else lv%sm2 => lv%sm end if - + if (.not.allocated(lv%aggr)) allocate(mld_c_base_aggregator_type :: lv%aggr,stat=info) + if (allocated(lv%aggr)) call lv%aggr%default() + return end subroutine c_base_onelev_default @@ -508,6 +528,14 @@ contains end if lvout%sm2 => lvout%sm end if + if (allocated(lv%aggr)) then + call lv%aggr%clone(lvout%aggr,info) + else + if (allocated(lvout%aggr)) then + call lvout%aggr%free(info) + if (info==psb_success_) deallocate(lvout%aggr,stat=info) + end if + end if if (info == psb_success_) call lv%parms%clone(lvout%parms,info) if (info == psb_success_) call lv%ac%clone(lvout%ac,info) if (info == psb_success_) call lv%tprol%clone(lvout%tprol,info) @@ -538,8 +566,9 @@ contains call move_alloc(lv%sm2a,b%sm2a) b%sm2 =>b%sm end if - - if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) + + call move_alloc(lv%aggr,b%aggr) + if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(lv%map,b%map,info) diff --git a/mlprec/mld_c_prec_mod.f90 b/mlprec/mld_c_prec_mod.f90 index 4ec6e9d2..a06680b8 100644 --- a/mlprec/mld_c_prec_mod.f90 +++ b/mlprec/mld_c_prec_mod.f90 @@ -55,7 +55,8 @@ module mld_c_prec_mod interface mld_precset module procedure mld_c_iprecsetsm, mld_c_iprecsetsv, & & mld_c_iprecseti, mld_c_iprecsetc, mld_c_iprecsetr, & - & mld_c_cprecseti, mld_c_cprecsetc, mld_c_cprecsetr + & mld_c_cprecseti, mld_c_cprecsetc, mld_c_cprecsetr, & + & mld_c_iprecsetag end interface mld_precset interface mld_extprol_bld @@ -97,6 +98,14 @@ contains call p%set(val,info, pos=pos) end subroutine mld_c_iprecsetsv + subroutine mld_c_iprecsetag(p,val,info,pos) + type(mld_cprec_type), intent(inout) :: p + class(mld_c_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + call p%set(val,info, pos=pos) + end subroutine mld_c_iprecsetag + subroutine mld_c_iprecseti(p,what,val,info,pos) type(mld_cprec_type), intent(inout) :: p integer(psb_ipk_), intent(in) :: what diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 index 7b0160d7..073f303c 100644 --- a/mlprec/mld_c_prec_type.f90 +++ b/mlprec/mld_c_prec_type.f90 @@ -55,6 +55,7 @@ module mld_c_prec_type use mld_base_prec_type use mld_c_base_solver_mod use mld_c_base_smoother_mod + use mld_c_base_aggregator_mod use mld_c_onelev_mod use psb_prec_mod, only : psb_cprec_type @@ -92,8 +93,8 @@ module mld_c_prec_type ! 2. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_levs = 20_psb_ipk_ ! 3. min_cr_ratio = 1.5 - real(psb_spk_) :: min_cr_ratio = 1.5_psb_spk_ - real(psb_spk_) :: op_complexity=szero + real(psb_spk_) :: min_cr_ratio = 1.5_psb_spk_ + real(psb_spk_) :: op_complexity = szero ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -126,6 +127,7 @@ module mld_c_prec_type procedure, pass(prec) :: sizeof => mld_cprec_sizeof procedure, pass(prec) :: setsm => mld_cprecsetsm procedure, pass(prec) :: setsv => mld_cprecsetsv + procedure, pass(prec) :: setag => mld_cprecsetag procedure, pass(prec) :: seti => mld_cprecseti procedure, pass(prec) :: setc => mld_cprecsetc procedure, pass(prec) :: setr => mld_cprecsetr @@ -133,7 +135,7 @@ module mld_c_prec_type procedure, pass(prec) :: csetc => mld_ccprecsetc procedure, pass(prec) :: csetr => mld_ccprecsetr generic, public :: set => seti, setc, setr, & - & cseti, csetc, csetr, setsm, setsv + & cseti, csetc, csetr, setsm, setsv, setag procedure, pass(prec) :: get_smoother => mld_c_get_smootherp procedure, pass(prec) :: get_solver => mld_c_get_solverp procedure, pass(prec) :: move_alloc => c_prec_move_alloc @@ -234,6 +236,15 @@ module mld_c_prec_type integer(psb_ipk_), optional, intent(in) :: ilev,ilmax character(len=*), optional, intent(in) :: pos end subroutine mld_cprecsetsv + subroutine mld_cprecsetag(prec,val,info,ilev,pos) + import :: psb_cspmat_type, psb_desc_type, psb_spk_, & + & mld_cprec_type, mld_c_base_aggregator_type, psb_ipk_ + class(mld_cprec_type), intent(inout) :: prec + class(mld_c_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev + character(len=*), optional, intent(in) :: pos + end subroutine mld_cprecsetag subroutine mld_cprecseti(prec,what,val,info,ilev,ilmax,pos) import :: psb_cspmat_type, psb_desc_type, psb_spk_, & & mld_cprec_type, psb_ipk_ diff --git a/mlprec/mld_c_symdec_aggregator_mod.f90 b/mlprec/mld_c_symdec_aggregator_mod.f90 new file mode 100644 index 00000000..4738c595 --- /dev/null +++ b/mlprec/mld_c_symdec_aggregator_mod.f90 @@ -0,0 +1,128 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This version only differs from the +! basic decoupled aggregation algorithm because it works on (the +! pattern of) A+A^T instead of A. +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_c_symdec_aggregator_mod + + use mld_c_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_c_base_aggregator_type) :: mld_c_symdec_aggregator_type + + contains + procedure, pass(ag) :: tprol => mld_c_symdec_aggregator_build_tprol + procedure, nopass :: fmt => mld_c_symdec_aggregator_fmt + end type mld_c_symdec_aggregator_type + + + interface + subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_c_symdec_aggregator_type, psb_desc_type, psb_cspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_c_symdec_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_symdec_aggregator_build_tprol + end interface + + +contains + + function mld_c_symdec_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Symmetric Decoupled aggregation" + end function mld_c_symdec_aggregator_fmt + +end module mld_c_symdec_aggregator_mod diff --git a/mlprec/mld_d_base_aggregator_mod.f90 b/mlprec/mld_d_base_aggregator_mod.f90 new file mode 100644 index 00000000..fd231c39 --- /dev/null +++ b/mlprec/mld_d_base_aggregator_mod.f90 @@ -0,0 +1,209 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. The basic version is the +! decoupled aggregation algorithm presented in +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_d_base_aggregator_mod + + use mld_base_prec_type, only : mld_dml_parms + use psb_base_mod, only : psb_dspmat_type, psb_d_vect_type, & + & psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & + & psb_erractionsave, psb_error_handler, psb_success_ + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type mld_d_base_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_d_base_aggregator_build_tprol + procedure, pass(ag) :: mat_asb => mld_d_base_aggregator_mat_asb + procedure, pass(ag) :: update_level => mld_d_base_aggregator_update_level + procedure, pass(ag) :: clone => mld_d_base_aggregator_clone + procedure, pass(ag) :: free => mld_d_base_aggregator_free + procedure, pass(ag) :: default => mld_d_base_aggregator_default + procedure, pass(ag) :: descr => mld_d_base_aggregator_descr + procedure, nopass :: fmt => mld_d_base_aggregator_fmt + end type mld_d_base_aggregator_type + + + interface + subroutine mld_d_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_d_base_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_d_base_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_base_aggregator_build_tprol + end interface + + interface + subroutine mld_d_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,& + & op_prol,op_restr,info) + import :: mld_d_base_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_d_base_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(psb_dspmat_type), intent(inout) :: op_prol + type(psb_dspmat_type), intent(out) :: ac,op_restr + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_base_aggregator_mat_asb + end interface + +contains + + subroutine mld_d_base_aggregator_update_level(ag,agnext,info) + implicit none + class(mld_d_base_aggregator_type), target, intent(inout) :: ag, agnext + integer(psb_ipk_), intent(out) :: info + + ! + ! Base version does nothing. + ! + info = 0 + end subroutine mld_d_base_aggregator_update_level + + subroutine mld_d_base_aggregator_clone(ag,agnext,info) + implicit none + class(mld_d_base_aggregator_type), intent(inout) :: ag + class(mld_d_base_aggregator_type), allocatable, intent(inout) :: agnext + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(agnext)) then + call agnext%free(info) + if (info == 0) deallocate(agnext,stat=info) + end if + if (info /= 0) return + allocate(agnext,source=ag,stat=info) + + end subroutine mld_d_base_aggregator_clone + + subroutine mld_d_base_aggregator_free(ag,info) + implicit none + class(mld_d_base_aggregator_type), intent(inout) :: ag + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + return + end subroutine mld_d_base_aggregator_free + + subroutine mld_d_base_aggregator_default(ag) + implicit none + class(mld_d_base_aggregator_type), intent(inout) :: ag + + ! Here we need do nothing + + return + end subroutine mld_d_base_aggregator_default + + function mld_d_base_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Decoupled aggregation" + end function mld_d_base_aggregator_fmt + + subroutine mld_d_base_aggregator_descr(ag,parms,iout,info) + implicit none + class(mld_d_base_aggregator_type), intent(in) :: ag + type(mld_dml_parms), intent(in) :: parms + integer(psb_ipk_), intent(in) :: iout + integer(psb_ipk_), intent(out) :: info + + call parms%mldescr(iout,info,aggr_name=ag%fmt()) + + return + end subroutine mld_d_base_aggregator_descr + +end module mld_d_base_aggregator_mod diff --git a/mlprec/mld_d_bcmatch_aggregator_mod.F90 b/mlprec/mld_d_bcmatch_aggregator_mod.F90 new file mode 100644 index 00000000..5f43ac91 --- /dev/null +++ b/mlprec/mld_d_bcmatch_aggregator_mod.F90 @@ -0,0 +1,262 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This variant is based on the hybrid method +! presented in +! +! S. Gratton, P. Henon, P. Jiranek and X. Vasseur: +! Reducing complexity of algebraic multigrid by aggregation +! Numerical Lin. Algebra with Applications, 2016, 23:501-518 +! + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + +module bcm_csr_type_mod + use iso_c_binding + type, bind(c):: bcm_Vector + type(c_ptr) :: data + integer(c_int) :: size + integer(c_int) :: owns_data + end type + + type, bind(c):: bcm_CSRMatrix + type(c_ptr) :: i + type(c_ptr) :: j + integer(c_int) :: num_rows + integer(c_int) :: num_cols + integer(c_int) :: num_nonzeros + integer(c_int) :: owns_data + type(c_ptr) :: data + end type +end module bcm_csr_type_mod + +module mld_d_bcmatch_aggregator_mod + use mld_d_base_aggregator_mod + use bcm_csr_type_mod + + type, extends(mld_d_base_aggregator_type) :: mld_d_bcmatch_aggregator_type + integer(psb_ipk_) :: matching_alg + integer(psb_ipk_) :: n_sweeps + real(psb_dpk_), allocatable :: w_tmp(:) + type(bcm_Vector) :: w_par + integer(psb_ipk_) :: max_csize + integer(psb_ipk_) :: max_nlevels + !type(psb_d_vect_type) :: w + contains + procedure, pass(ag) :: bld_tprol => mld_d_bcmatch_aggregator_build_tprol + procedure, pass(ag) :: set => d_bcmatch_aggr_cseti + procedure, pass(ag) :: default =>d_bcmatch_aggr_set_default + procedure, pass(ag) :: mat_asb => mld_d_bcmatch_aggregator_mat_asb + procedure, pass(ag) :: update_level => d_bcmatch_aggregator_update_level +!!$ procedure, pass(ag) :: clone => mld_d_base_aggregator_clone +!!$ procedure, pass(ag) :: free => mld_d_bcmatch_aggregator_free +!!$ procedure, pass(ag) :: default => mld_d_base_aggregator_default + procedure, nopass :: fmt => mld_d_bcmatch_aggregator_fmt + end type mld_d_bcmatch_aggregator_type + + + interface + subroutine mld_d_bcmatch_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_d_bcmatch_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_bcmatch_aggregator_build_tprol + end interface + + interface + subroutine mld_d_bcmatch_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,& + & op_prol,op_restr,info) + import :: mld_d_bcmatch_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(psb_dspmat_type), intent(inout) :: op_prol + type(psb_dspmat_type), intent(out) :: ac,op_restr + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_bcmatch_aggregator_mat_asb + end interface + + +contains + + + function mld_d_bcmatch_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "BootCMatch aggregation" + end function mld_d_bcmatch_aggregator_fmt + + subroutine d_bcmatch_aggregator_update_level(ag,agnext,info) + implicit none + class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag + class(mld_d_base_aggregator_type), target, intent(inout) :: agnext + integer(psb_ipk_), intent(out) :: info + + ! + ! + select type(agnext) + type is (mld_d_bcmatch_aggregator_type) + agnext%matching_alg=ag%matching_alg + agnext%n_sweeps=ag%n_sweeps + agnext%max_csize=ag%max_csize + agnext%max_nlevels=ag%max_nlevels + agnext%w_par=ag%w_par + end select + info = 0 + end subroutine d_bcmatch_aggregator_update_level + + subroutine d_bcmatch_aggr_cseti(ag,what,val,info) + + Implicit None + + ! Arguments + class(mld_d_bcmatch_aggregator_type), intent(inout) :: ag + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, iwhat + character(len=20) :: name='d_bcmatch_aggr_cseti' + info = psb_success_ + + select case(what) + case('BCM_MATCH_ALG') + ag%matching_alg=val + case('BCM_SWEEPS') + ag%n_sweeps=val + case('BCM_MAX_CSIZE') + ag%max_csize=val + case('BCM_MAX_NLEVELS') + ag%max_nlevels=val + case('BCM_W_SIZE') + ag%w_par%size=val + ag%w_par%owns_data=0 + allocate(ag%w_tmp(val)) + ag%w_tmp = 1.0_psb_dpk_ + call set_cloc(ag%w_tmp, ag%w_par) + case default + end select + return + contains + subroutine set_cloc(vect,w_par) + real(psb_dpk_), target :: vect(:) + type(bcm_Vector) :: w_par + + w_par%data = c_loc(vect) + end subroutine set_cloc + + end subroutine d_bcmatch_aggr_cseti + + subroutine d_bcmatch_aggr_set_default(ag) + + Implicit None + + ! Arguments + class(mld_d_bcmatch_aggregator_type), intent(inout) :: ag + character(len=20) :: name='d_bcmatch_aggr_set_default' + ag%matching_alg=0 + ag%n_sweeps=1 + ag%max_nlevels=36 + ag%max_csize=10 + + return + + end subroutine d_bcmatch_aggr_set_default + +!!$ subroutine d_bcmatch_aggregator_free(ag,info) +!!$ implicit none +!!$ class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ info = 0 +!!$ end subroutine d_bcmatch_aggregator_free + + +end module mld_d_bcmatch_aggregator_mod diff --git a/mlprec/mld_d_hybrid_aggregator_mod.F90 b/mlprec/mld_d_hybrid_aggregator_mod.F90 new file mode 100644 index 00000000..50e3215c --- /dev/null +++ b/mlprec/mld_d_hybrid_aggregator_mod.F90 @@ -0,0 +1,131 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This variant is based on the hybrid method +! presented in +! +! S. Gratton, P. Henon, P. Jiranek and X. Vasseur: +! Reducing complexity of algebraic multigrid by aggregation +! Numerical Lin. Algebra with Applications, 2016, 23:501-518 +! +module mld_d_hybrid_aggregator_mod + + use mld_d_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_d_base_aggregator_type) :: mld_d_hybrid_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_d_hybrid_aggregator_build_tprol +!!$ procedure, pass(ag) :: mat_asb => mld_d_base_aggregator_mat_asb +!!$ procedure, pass(ag) :: update_level => mld_d_base_aggregator_update_level +!!$ procedure, pass(ag) :: clone => mld_d_base_aggregator_clone +!!$ procedure, pass(ag) :: free => mld_d_base_aggregator_free +!!$ procedure, pass(ag) :: default => mld_d_base_aggregator_default + procedure, nopass :: fmt => mld_d_hybrid_aggregator_fmt + end type mld_d_hybrid_aggregator_type + + + interface + subroutine mld_d_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_d_hybrid_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_d_hybrid_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_hybrid_aggregator_build_tprol + end interface + +contains + + + function mld_d_hybrid_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Hybrid Decoupled aggregation" + end function mld_d_hybrid_aggregator_fmt + + +end module mld_d_hybrid_aggregator_mod diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index e73856a0..810d5831 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -137,6 +137,31 @@ module mld_d_inner_mod end subroutine mld_d_dec_map_bld end interface mld_dec_map_bld + interface mld_hyb_map_bld + subroutine mld_d_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: theta + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_hyb_map_bld + end interface mld_hyb_map_bld + + interface mld_map_to_tprol + subroutine mld_d_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + use mld_d_prec_type, only : mld_d_onelev_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_map_to_tprol + end interface mld_map_to_tprol + interface mld_lev_mat_asb subroutine mld_d_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index 4420f870..30c7376c 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -55,6 +55,7 @@ module mld_d_onelev_mod use mld_base_prec_type use mld_d_base_smoother_mod + use mld_d_base_aggregator_mod use psb_base_mod, only : psb_dspmat_type, psb_d_vect_type, & & psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, & & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & @@ -136,9 +137,10 @@ module mld_d_onelev_mod & d_wrk_clone, d_wrk_move_alloc, d_wrk_cnv type mld_d_onelev_type - class(mld_d_base_smoother_type), allocatable :: sm, sm2a - class(mld_d_base_smoother_type), pointer :: sm2 => null() - class(mld_dmlprec_wrk_type), allocatable :: wrk + class(mld_d_base_smoother_type), allocatable :: sm, sm2a + class(mld_d_base_smoother_type), pointer :: sm2 => null() + class(mld_dmlprec_wrk_type), allocatable :: wrk + class(mld_d_base_aggregator_type), allocatable :: aggr type(mld_dml_parms) :: parms type(psb_dspmat_type) :: ac integer(psb_ipk_) :: ac_nz_loc, ac_nz_tot @@ -166,8 +168,9 @@ module mld_d_onelev_mod procedure, pass(lv) :: csetc => mld_d_base_onelev_csetc procedure, pass(lv) :: setsm => mld_d_base_onelev_setsm procedure, pass(lv) :: setsv => mld_d_base_onelev_setsv + procedure, pass(lv) :: setag => mld_d_base_onelev_setag generic, public :: set => seti, setr, setc, & - & cseti, csetr, csetc, setsm, setsv + & cseti, csetr, csetc, setsm, setsv, setag procedure, pass(lv) :: sizeof => d_base_onelev_sizeof procedure, pass(lv) :: get_nzeros => d_base_onelev_get_nzeros procedure, pass(lv) :: get_wrksz => d_base_onelev_get_wrksize @@ -299,6 +302,20 @@ module mld_d_onelev_mod end subroutine mld_d_base_onelev_setsv end interface + interface + subroutine mld_d_base_onelev_setag(lv,val,info,pos) + import :: psb_dpk_, mld_d_onelev_type, mld_d_base_aggregator_type, & + & psb_ipk_, psb_long_int_k_, psb_desc_type + Implicit None + + ! Arguments + class(mld_d_onelev_type), target, intent(inout) :: lv + class(mld_d_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + end subroutine mld_d_base_onelev_setag + end interface + interface subroutine mld_d_base_onelev_setc(lv,what,val,info,pos) import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & @@ -323,9 +340,9 @@ module mld_d_onelev_mod class(mld_d_onelev_type), intent(inout) :: lv integer(psb_ipk_), intent(in) :: what - real(psb_dpk_), intent(in) :: val + real(psb_dpk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info - character(len=*), optional, intent(in) :: pos + character(len=*), optional, intent(in) :: pos end subroutine mld_d_base_onelev_setr end interface @@ -450,7 +467,8 @@ contains Implicit None ! Arguments - class(mld_d_onelev_type), target, intent(inout) :: lv + class(mld_d_onelev_type), target, intent(inout) :: lv + integer(psb_ipk_) :: info lv%parms%sweeps_pre = 1 lv%parms%sweeps_post = 1 @@ -473,7 +491,9 @@ contains else lv%sm2 => lv%sm end if - + if (.not.allocated(lv%aggr)) allocate(mld_d_base_aggregator_type :: lv%aggr,stat=info) + if (allocated(lv%aggr)) call lv%aggr%default() + return end subroutine d_base_onelev_default @@ -508,6 +528,14 @@ contains end if lvout%sm2 => lvout%sm end if + if (allocated(lv%aggr)) then + call lv%aggr%clone(lvout%aggr,info) + else + if (allocated(lvout%aggr)) then + call lvout%aggr%free(info) + if (info==psb_success_) deallocate(lvout%aggr,stat=info) + end if + end if if (info == psb_success_) call lv%parms%clone(lvout%parms,info) if (info == psb_success_) call lv%ac%clone(lvout%ac,info) if (info == psb_success_) call lv%tprol%clone(lvout%tprol,info) @@ -538,8 +566,9 @@ contains call move_alloc(lv%sm2a,b%sm2a) b%sm2 =>b%sm end if - - if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) + + call move_alloc(lv%aggr,b%aggr) + if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(lv%map,b%map,info) diff --git a/mlprec/mld_d_prec_mod.f90 b/mlprec/mld_d_prec_mod.f90 index 152c860f..7e695111 100644 --- a/mlprec/mld_d_prec_mod.f90 +++ b/mlprec/mld_d_prec_mod.f90 @@ -55,7 +55,8 @@ module mld_d_prec_mod interface mld_precset module procedure mld_d_iprecsetsm, mld_d_iprecsetsv, & & mld_d_iprecseti, mld_d_iprecsetc, mld_d_iprecsetr, & - & mld_d_cprecseti, mld_d_cprecsetc, mld_d_cprecsetr + & mld_d_cprecseti, mld_d_cprecsetc, mld_d_cprecsetr, & + & mld_d_iprecsetag end interface mld_precset interface mld_extprol_bld @@ -97,6 +98,14 @@ contains call p%set(val,info, pos=pos) end subroutine mld_d_iprecsetsv + subroutine mld_d_iprecsetag(p,val,info,pos) + type(mld_dprec_type), intent(inout) :: p + class(mld_d_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + call p%set(val,info, pos=pos) + end subroutine mld_d_iprecsetag + subroutine mld_d_iprecseti(p,what,val,info,pos) type(mld_dprec_type), intent(inout) :: p integer(psb_ipk_), intent(in) :: what diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index 1aa9409c..f62f4760 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -55,6 +55,7 @@ module mld_d_prec_type use mld_base_prec_type use mld_d_base_solver_mod use mld_d_base_smoother_mod + use mld_d_base_aggregator_mod use mld_d_onelev_mod use psb_prec_mod, only : psb_dprec_type @@ -92,8 +93,8 @@ module mld_d_prec_type ! 2. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_levs = 20_psb_ipk_ ! 3. min_cr_ratio = 1.5 - real(psb_dpk_) :: min_cr_ratio = 1.5_psb_dpk_ - real(psb_dpk_) :: op_complexity=dzero + real(psb_dpk_) :: min_cr_ratio = 1.5_psb_dpk_ + real(psb_dpk_) :: op_complexity = dzero ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -126,6 +127,7 @@ module mld_d_prec_type procedure, pass(prec) :: sizeof => mld_dprec_sizeof procedure, pass(prec) :: setsm => mld_dprecsetsm procedure, pass(prec) :: setsv => mld_dprecsetsv + procedure, pass(prec) :: setag => mld_dprecsetag procedure, pass(prec) :: seti => mld_dprecseti procedure, pass(prec) :: setc => mld_dprecsetc procedure, pass(prec) :: setr => mld_dprecsetr @@ -133,7 +135,7 @@ module mld_d_prec_type procedure, pass(prec) :: csetc => mld_dcprecsetc procedure, pass(prec) :: csetr => mld_dcprecsetr generic, public :: set => seti, setc, setr, & - & cseti, csetc, csetr, setsm, setsv + & cseti, csetc, csetr, setsm, setsv, setag procedure, pass(prec) :: get_smoother => mld_d_get_smootherp procedure, pass(prec) :: get_solver => mld_d_get_solverp procedure, pass(prec) :: move_alloc => d_prec_move_alloc @@ -234,6 +236,15 @@ module mld_d_prec_type integer(psb_ipk_), optional, intent(in) :: ilev,ilmax character(len=*), optional, intent(in) :: pos end subroutine mld_dprecsetsv + subroutine mld_dprecsetag(prec,val,info,ilev,pos) + import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & + & mld_dprec_type, mld_d_base_aggregator_type, psb_ipk_ + class(mld_dprec_type), intent(inout) :: prec + class(mld_d_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev + character(len=*), optional, intent(in) :: pos + end subroutine mld_dprecsetag subroutine mld_dprecseti(prec,what,val,info,ilev,ilmax,pos) import :: psb_dspmat_type, psb_desc_type, psb_dpk_, & & mld_dprec_type, psb_ipk_ diff --git a/mlprec/mld_d_symdec_aggregator_mod.f90 b/mlprec/mld_d_symdec_aggregator_mod.f90 new file mode 100644 index 00000000..2532fac9 --- /dev/null +++ b/mlprec/mld_d_symdec_aggregator_mod.f90 @@ -0,0 +1,128 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This version only differs from the +! basic decoupled aggregation algorithm because it works on (the +! pattern of) A+A^T instead of A. +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_d_symdec_aggregator_mod + + use mld_d_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_d_base_aggregator_type) :: mld_d_symdec_aggregator_type + + contains + procedure, pass(ag) :: tprol => mld_d_symdec_aggregator_build_tprol + procedure, nopass :: fmt => mld_d_symdec_aggregator_fmt + end type mld_d_symdec_aggregator_type + + + interface + subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_d_symdec_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_d_symdec_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_symdec_aggregator_build_tprol + end interface + + +contains + + function mld_d_symdec_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Symmetric Decoupled aggregation" + end function mld_d_symdec_aggregator_fmt + +end module mld_d_symdec_aggregator_mod diff --git a/mlprec/mld_s_base_aggregator_mod.f90 b/mlprec/mld_s_base_aggregator_mod.f90 new file mode 100644 index 00000000..3021e19c --- /dev/null +++ b/mlprec/mld_s_base_aggregator_mod.f90 @@ -0,0 +1,209 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. The basic version is the +! decoupled aggregation algorithm presented in +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_s_base_aggregator_mod + + use mld_base_prec_type, only : mld_sml_parms + use psb_base_mod, only : psb_sspmat_type, psb_s_vect_type, & + & psb_s_base_vect_type, psb_slinmap_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & + & psb_erractionsave, psb_error_handler, psb_success_ + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type mld_s_base_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_s_base_aggregator_build_tprol + procedure, pass(ag) :: mat_asb => mld_s_base_aggregator_mat_asb + procedure, pass(ag) :: update_level => mld_s_base_aggregator_update_level + procedure, pass(ag) :: clone => mld_s_base_aggregator_clone + procedure, pass(ag) :: free => mld_s_base_aggregator_free + procedure, pass(ag) :: default => mld_s_base_aggregator_default + procedure, pass(ag) :: descr => mld_s_base_aggregator_descr + procedure, nopass :: fmt => mld_s_base_aggregator_fmt + end type mld_s_base_aggregator_type + + + interface + subroutine mld_s_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_s_base_aggregator_type, psb_desc_type, psb_sspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_s_base_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_base_aggregator_build_tprol + end interface + + interface + subroutine mld_s_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,& + & op_prol,op_restr,info) + import :: mld_s_base_aggregator_type, psb_desc_type, psb_sspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_s_base_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(psb_sspmat_type), intent(inout) :: op_prol + type(psb_sspmat_type), intent(out) :: ac,op_restr + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_base_aggregator_mat_asb + end interface + +contains + + subroutine mld_s_base_aggregator_update_level(ag,agnext,info) + implicit none + class(mld_s_base_aggregator_type), target, intent(inout) :: ag, agnext + integer(psb_ipk_), intent(out) :: info + + ! + ! Base version does nothing. + ! + info = 0 + end subroutine mld_s_base_aggregator_update_level + + subroutine mld_s_base_aggregator_clone(ag,agnext,info) + implicit none + class(mld_s_base_aggregator_type), intent(inout) :: ag + class(mld_s_base_aggregator_type), allocatable, intent(inout) :: agnext + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(agnext)) then + call agnext%free(info) + if (info == 0) deallocate(agnext,stat=info) + end if + if (info /= 0) return + allocate(agnext,source=ag,stat=info) + + end subroutine mld_s_base_aggregator_clone + + subroutine mld_s_base_aggregator_free(ag,info) + implicit none + class(mld_s_base_aggregator_type), intent(inout) :: ag + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + return + end subroutine mld_s_base_aggregator_free + + subroutine mld_s_base_aggregator_default(ag) + implicit none + class(mld_s_base_aggregator_type), intent(inout) :: ag + + ! Here we need do nothing + + return + end subroutine mld_s_base_aggregator_default + + function mld_s_base_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Decoupled aggregation" + end function mld_s_base_aggregator_fmt + + subroutine mld_s_base_aggregator_descr(ag,parms,iout,info) + implicit none + class(mld_s_base_aggregator_type), intent(in) :: ag + type(mld_sml_parms), intent(in) :: parms + integer(psb_ipk_), intent(in) :: iout + integer(psb_ipk_), intent(out) :: info + + call parms%mldescr(iout,info,aggr_name=ag%fmt()) + + return + end subroutine mld_s_base_aggregator_descr + +end module mld_s_base_aggregator_mod diff --git a/mlprec/mld_s_hybrid_aggregator_mod.F90 b/mlprec/mld_s_hybrid_aggregator_mod.F90 new file mode 100644 index 00000000..5db38098 --- /dev/null +++ b/mlprec/mld_s_hybrid_aggregator_mod.F90 @@ -0,0 +1,131 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This variant is based on the hybrid method +! presented in +! +! S. Gratton, P. Henon, P. Jiranek and X. Vasseur: +! Reducing complexity of algebraic multigrid by aggregation +! Numerical Lin. Algebra with Applications, 2016, 23:501-518 +! +module mld_s_hybrid_aggregator_mod + + use mld_s_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_s_base_aggregator_type) :: mld_s_hybrid_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_s_hybrid_aggregator_build_tprol +!!$ procedure, pass(ag) :: mat_asb => mld_s_base_aggregator_mat_asb +!!$ procedure, pass(ag) :: update_level => mld_s_base_aggregator_update_level +!!$ procedure, pass(ag) :: clone => mld_s_base_aggregator_clone +!!$ procedure, pass(ag) :: free => mld_s_base_aggregator_free +!!$ procedure, pass(ag) :: default => mld_s_base_aggregator_default + procedure, nopass :: fmt => mld_s_hybrid_aggregator_fmt + end type mld_s_hybrid_aggregator_type + + + interface + subroutine mld_s_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_s_hybrid_aggregator_type, psb_desc_type, psb_sspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_s_hybrid_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_hybrid_aggregator_build_tprol + end interface + +contains + + + function mld_s_hybrid_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Hybrid Decoupled aggregation" + end function mld_s_hybrid_aggregator_fmt + + +end module mld_s_hybrid_aggregator_mod diff --git a/mlprec/mld_s_inner_mod.f90 b/mlprec/mld_s_inner_mod.f90 index 581b943e..73baa6a8 100644 --- a/mlprec/mld_s_inner_mod.f90 +++ b/mlprec/mld_s_inner_mod.f90 @@ -137,6 +137,31 @@ module mld_s_inner_mod end subroutine mld_s_dec_map_bld end interface mld_dec_map_bld + interface mld_hyb_map_bld + subroutine mld_s_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: theta + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_hyb_map_bld + end interface mld_hyb_map_bld + + interface mld_map_to_tprol + subroutine mld_s_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + use mld_s_prec_type, only : mld_s_onelev_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_map_to_tprol + end interface mld_map_to_tprol + interface mld_lev_mat_asb subroutine mld_s_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) diff --git a/mlprec/mld_s_onelev_mod.f90 b/mlprec/mld_s_onelev_mod.f90 index e135b0f4..baba63cc 100644 --- a/mlprec/mld_s_onelev_mod.f90 +++ b/mlprec/mld_s_onelev_mod.f90 @@ -55,6 +55,7 @@ module mld_s_onelev_mod use mld_base_prec_type use mld_s_base_smoother_mod + use mld_s_base_aggregator_mod use psb_base_mod, only : psb_sspmat_type, psb_s_vect_type, & & psb_s_base_vect_type, psb_slinmap_type, psb_spk_, & & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & @@ -136,9 +137,10 @@ module mld_s_onelev_mod & s_wrk_clone, s_wrk_move_alloc, s_wrk_cnv type mld_s_onelev_type - class(mld_s_base_smoother_type), allocatable :: sm, sm2a - class(mld_s_base_smoother_type), pointer :: sm2 => null() - class(mld_smlprec_wrk_type), allocatable :: wrk + class(mld_s_base_smoother_type), allocatable :: sm, sm2a + class(mld_s_base_smoother_type), pointer :: sm2 => null() + class(mld_smlprec_wrk_type), allocatable :: wrk + class(mld_s_base_aggregator_type), allocatable :: aggr type(mld_sml_parms) :: parms type(psb_sspmat_type) :: ac integer(psb_ipk_) :: ac_nz_loc, ac_nz_tot @@ -166,8 +168,9 @@ module mld_s_onelev_mod procedure, pass(lv) :: csetc => mld_s_base_onelev_csetc procedure, pass(lv) :: setsm => mld_s_base_onelev_setsm procedure, pass(lv) :: setsv => mld_s_base_onelev_setsv + procedure, pass(lv) :: setag => mld_s_base_onelev_setag generic, public :: set => seti, setr, setc, & - & cseti, csetr, csetc, setsm, setsv + & cseti, csetr, csetc, setsm, setsv, setag procedure, pass(lv) :: sizeof => s_base_onelev_sizeof procedure, pass(lv) :: get_nzeros => s_base_onelev_get_nzeros procedure, pass(lv) :: get_wrksz => s_base_onelev_get_wrksize @@ -299,6 +302,20 @@ module mld_s_onelev_mod end subroutine mld_s_base_onelev_setsv end interface + interface + subroutine mld_s_base_onelev_setag(lv,val,info,pos) + import :: psb_spk_, mld_s_onelev_type, mld_s_base_aggregator_type, & + & psb_ipk_, psb_long_int_k_, psb_desc_type + Implicit None + + ! Arguments + class(mld_s_onelev_type), target, intent(inout) :: lv + class(mld_s_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + end subroutine mld_s_base_onelev_setag + end interface + interface subroutine mld_s_base_onelev_setc(lv,what,val,info,pos) import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & @@ -323,9 +340,9 @@ module mld_s_onelev_mod class(mld_s_onelev_type), intent(inout) :: lv integer(psb_ipk_), intent(in) :: what - real(psb_spk_), intent(in) :: val + real(psb_spk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info - character(len=*), optional, intent(in) :: pos + character(len=*), optional, intent(in) :: pos end subroutine mld_s_base_onelev_setr end interface @@ -450,7 +467,8 @@ contains Implicit None ! Arguments - class(mld_s_onelev_type), target, intent(inout) :: lv + class(mld_s_onelev_type), target, intent(inout) :: lv + integer(psb_ipk_) :: info lv%parms%sweeps_pre = 1 lv%parms%sweeps_post = 1 @@ -473,7 +491,9 @@ contains else lv%sm2 => lv%sm end if - + if (.not.allocated(lv%aggr)) allocate(mld_s_base_aggregator_type :: lv%aggr,stat=info) + if (allocated(lv%aggr)) call lv%aggr%default() + return end subroutine s_base_onelev_default @@ -508,6 +528,14 @@ contains end if lvout%sm2 => lvout%sm end if + if (allocated(lv%aggr)) then + call lv%aggr%clone(lvout%aggr,info) + else + if (allocated(lvout%aggr)) then + call lvout%aggr%free(info) + if (info==psb_success_) deallocate(lvout%aggr,stat=info) + end if + end if if (info == psb_success_) call lv%parms%clone(lvout%parms,info) if (info == psb_success_) call lv%ac%clone(lvout%ac,info) if (info == psb_success_) call lv%tprol%clone(lvout%tprol,info) @@ -538,8 +566,9 @@ contains call move_alloc(lv%sm2a,b%sm2a) b%sm2 =>b%sm end if - - if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) + + call move_alloc(lv%aggr,b%aggr) + if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(lv%map,b%map,info) diff --git a/mlprec/mld_s_prec_mod.f90 b/mlprec/mld_s_prec_mod.f90 index b7e35d27..27e9c75a 100644 --- a/mlprec/mld_s_prec_mod.f90 +++ b/mlprec/mld_s_prec_mod.f90 @@ -55,7 +55,8 @@ module mld_s_prec_mod interface mld_precset module procedure mld_s_iprecsetsm, mld_s_iprecsetsv, & & mld_s_iprecseti, mld_s_iprecsetc, mld_s_iprecsetr, & - & mld_s_cprecseti, mld_s_cprecsetc, mld_s_cprecsetr + & mld_s_cprecseti, mld_s_cprecsetc, mld_s_cprecsetr, & + & mld_s_iprecsetag end interface mld_precset interface mld_extprol_bld @@ -97,6 +98,14 @@ contains call p%set(val,info, pos=pos) end subroutine mld_s_iprecsetsv + subroutine mld_s_iprecsetag(p,val,info,pos) + type(mld_sprec_type), intent(inout) :: p + class(mld_s_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + call p%set(val,info, pos=pos) + end subroutine mld_s_iprecsetag + subroutine mld_s_iprecseti(p,what,val,info,pos) type(mld_sprec_type), intent(inout) :: p integer(psb_ipk_), intent(in) :: what diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 index 7b67ad66..a53c65e8 100644 --- a/mlprec/mld_s_prec_type.f90 +++ b/mlprec/mld_s_prec_type.f90 @@ -55,6 +55,7 @@ module mld_s_prec_type use mld_base_prec_type use mld_s_base_solver_mod use mld_s_base_smoother_mod + use mld_s_base_aggregator_mod use mld_s_onelev_mod use psb_prec_mod, only : psb_sprec_type @@ -92,8 +93,8 @@ module mld_s_prec_type ! 2. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_levs = 20_psb_ipk_ ! 3. min_cr_ratio = 1.5 - real(psb_spk_) :: min_cr_ratio = 1.5_psb_spk_ - real(psb_spk_) :: op_complexity=szero + real(psb_spk_) :: min_cr_ratio = 1.5_psb_spk_ + real(psb_spk_) :: op_complexity = szero ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -126,6 +127,7 @@ module mld_s_prec_type procedure, pass(prec) :: sizeof => mld_sprec_sizeof procedure, pass(prec) :: setsm => mld_sprecsetsm procedure, pass(prec) :: setsv => mld_sprecsetsv + procedure, pass(prec) :: setag => mld_sprecsetag procedure, pass(prec) :: seti => mld_sprecseti procedure, pass(prec) :: setc => mld_sprecsetc procedure, pass(prec) :: setr => mld_sprecsetr @@ -133,7 +135,7 @@ module mld_s_prec_type procedure, pass(prec) :: csetc => mld_scprecsetc procedure, pass(prec) :: csetr => mld_scprecsetr generic, public :: set => seti, setc, setr, & - & cseti, csetc, csetr, setsm, setsv + & cseti, csetc, csetr, setsm, setsv, setag procedure, pass(prec) :: get_smoother => mld_s_get_smootherp procedure, pass(prec) :: get_solver => mld_s_get_solverp procedure, pass(prec) :: move_alloc => s_prec_move_alloc @@ -234,6 +236,15 @@ module mld_s_prec_type integer(psb_ipk_), optional, intent(in) :: ilev,ilmax character(len=*), optional, intent(in) :: pos end subroutine mld_sprecsetsv + subroutine mld_sprecsetag(prec,val,info,ilev,pos) + import :: psb_sspmat_type, psb_desc_type, psb_spk_, & + & mld_sprec_type, mld_s_base_aggregator_type, psb_ipk_ + class(mld_sprec_type), intent(inout) :: prec + class(mld_s_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev + character(len=*), optional, intent(in) :: pos + end subroutine mld_sprecsetag subroutine mld_sprecseti(prec,what,val,info,ilev,ilmax,pos) import :: psb_sspmat_type, psb_desc_type, psb_spk_, & & mld_sprec_type, psb_ipk_ diff --git a/mlprec/mld_s_symdec_aggregator_mod.f90 b/mlprec/mld_s_symdec_aggregator_mod.f90 new file mode 100644 index 00000000..fbc89303 --- /dev/null +++ b/mlprec/mld_s_symdec_aggregator_mod.f90 @@ -0,0 +1,128 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This version only differs from the +! basic decoupled aggregation algorithm because it works on (the +! pattern of) A+A^T instead of A. +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_s_symdec_aggregator_mod + + use mld_s_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_s_base_aggregator_type) :: mld_s_symdec_aggregator_type + + contains + procedure, pass(ag) :: tprol => mld_s_symdec_aggregator_build_tprol + procedure, nopass :: fmt => mld_s_symdec_aggregator_fmt + end type mld_s_symdec_aggregator_type + + + interface + subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_s_symdec_aggregator_type, psb_desc_type, psb_sspmat_type, psb_spk_, & + & psb_ipk_, psb_long_int_k_, mld_sml_parms + implicit none + class(mld_s_symdec_aggregator_type), target, intent(inout) :: ag + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_symdec_aggregator_build_tprol + end interface + + +contains + + function mld_s_symdec_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Symmetric Decoupled aggregation" + end function mld_s_symdec_aggregator_fmt + +end module mld_s_symdec_aggregator_mod diff --git a/mlprec/mld_z_base_aggregator_mod.f90 b/mlprec/mld_z_base_aggregator_mod.f90 new file mode 100644 index 00000000..6f3a8096 --- /dev/null +++ b/mlprec/mld_z_base_aggregator_mod.f90 @@ -0,0 +1,209 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. The basic version is the +! decoupled aggregation algorithm presented in +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_z_base_aggregator_mod + + use mld_base_prec_type, only : mld_dml_parms + use psb_base_mod, only : psb_zspmat_type, psb_z_vect_type, & + & psb_z_base_vect_type, psb_zlinmap_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & + & psb_erractionsave, psb_error_handler, psb_success_ + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type mld_z_base_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_z_base_aggregator_build_tprol + procedure, pass(ag) :: mat_asb => mld_z_base_aggregator_mat_asb + procedure, pass(ag) :: update_level => mld_z_base_aggregator_update_level + procedure, pass(ag) :: clone => mld_z_base_aggregator_clone + procedure, pass(ag) :: free => mld_z_base_aggregator_free + procedure, pass(ag) :: default => mld_z_base_aggregator_default + procedure, pass(ag) :: descr => mld_z_base_aggregator_descr + procedure, nopass :: fmt => mld_z_base_aggregator_fmt + end type mld_z_base_aggregator_type + + + interface + subroutine mld_z_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_z_base_aggregator_type, psb_desc_type, psb_zspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_z_base_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_base_aggregator_build_tprol + end interface + + interface + subroutine mld_z_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,& + & op_prol,op_restr,info) + import :: mld_z_base_aggregator_type, psb_desc_type, psb_zspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_z_base_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(psb_zspmat_type), intent(inout) :: op_prol + type(psb_zspmat_type), intent(out) :: ac,op_restr + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_base_aggregator_mat_asb + end interface + +contains + + subroutine mld_z_base_aggregator_update_level(ag,agnext,info) + implicit none + class(mld_z_base_aggregator_type), target, intent(inout) :: ag, agnext + integer(psb_ipk_), intent(out) :: info + + ! + ! Base version does nothing. + ! + info = 0 + end subroutine mld_z_base_aggregator_update_level + + subroutine mld_z_base_aggregator_clone(ag,agnext,info) + implicit none + class(mld_z_base_aggregator_type), intent(inout) :: ag + class(mld_z_base_aggregator_type), allocatable, intent(inout) :: agnext + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(agnext)) then + call agnext%free(info) + if (info == 0) deallocate(agnext,stat=info) + end if + if (info /= 0) return + allocate(agnext,source=ag,stat=info) + + end subroutine mld_z_base_aggregator_clone + + subroutine mld_z_base_aggregator_free(ag,info) + implicit none + class(mld_z_base_aggregator_type), intent(inout) :: ag + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + return + end subroutine mld_z_base_aggregator_free + + subroutine mld_z_base_aggregator_default(ag) + implicit none + class(mld_z_base_aggregator_type), intent(inout) :: ag + + ! Here we need do nothing + + return + end subroutine mld_z_base_aggregator_default + + function mld_z_base_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Decoupled aggregation" + end function mld_z_base_aggregator_fmt + + subroutine mld_z_base_aggregator_descr(ag,parms,iout,info) + implicit none + class(mld_z_base_aggregator_type), intent(in) :: ag + type(mld_dml_parms), intent(in) :: parms + integer(psb_ipk_), intent(in) :: iout + integer(psb_ipk_), intent(out) :: info + + call parms%mldescr(iout,info,aggr_name=ag%fmt()) + + return + end subroutine mld_z_base_aggregator_descr + +end module mld_z_base_aggregator_mod diff --git a/mlprec/mld_z_hybrid_aggregator_mod.F90 b/mlprec/mld_z_hybrid_aggregator_mod.F90 new file mode 100644 index 00000000..be353fc4 --- /dev/null +++ b/mlprec/mld_z_hybrid_aggregator_mod.F90 @@ -0,0 +1,131 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This variant is based on the hybrid method +! presented in +! +! S. Gratton, P. Henon, P. Jiranek and X. Vasseur: +! Reducing complexity of algebraic multigrid by aggregation +! Numerical Lin. Algebra with Applications, 2016, 23:501-518 +! +module mld_z_hybrid_aggregator_mod + + use mld_z_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_z_base_aggregator_type) :: mld_z_hybrid_aggregator_type + + contains + procedure, pass(ag) :: bld_tprol => mld_z_hybrid_aggregator_build_tprol +!!$ procedure, pass(ag) :: mat_asb => mld_z_base_aggregator_mat_asb +!!$ procedure, pass(ag) :: update_level => mld_z_base_aggregator_update_level +!!$ procedure, pass(ag) :: clone => mld_z_base_aggregator_clone +!!$ procedure, pass(ag) :: free => mld_z_base_aggregator_free +!!$ procedure, pass(ag) :: default => mld_z_base_aggregator_default + procedure, nopass :: fmt => mld_z_hybrid_aggregator_fmt + end type mld_z_hybrid_aggregator_type + + + interface + subroutine mld_z_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_z_hybrid_aggregator_type, psb_desc_type, psb_zspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_z_hybrid_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_hybrid_aggregator_build_tprol + end interface + +contains + + + function mld_z_hybrid_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Hybrid Decoupled aggregation" + end function mld_z_hybrid_aggregator_fmt + + +end module mld_z_hybrid_aggregator_mod diff --git a/mlprec/mld_z_inner_mod.f90 b/mlprec/mld_z_inner_mod.f90 index 12860b7b..e6dcf6ee 100644 --- a/mlprec/mld_z_inner_mod.f90 +++ b/mlprec/mld_z_inner_mod.f90 @@ -137,6 +137,31 @@ module mld_z_inner_mod end subroutine mld_z_dec_map_bld end interface mld_dec_map_bld + interface mld_hyb_map_bld + subroutine mld_z_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: theta + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_hyb_map_bld + end interface mld_hyb_map_bld + + interface mld_map_to_tprol + subroutine mld_z_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + use mld_z_prec_type, only : mld_z_onelev_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_map_to_tprol + end interface mld_map_to_tprol + interface mld_lev_mat_asb subroutine mld_z_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) diff --git a/mlprec/mld_z_onelev_mod.f90 b/mlprec/mld_z_onelev_mod.f90 index 986ea076..be8eaaf3 100644 --- a/mlprec/mld_z_onelev_mod.f90 +++ b/mlprec/mld_z_onelev_mod.f90 @@ -55,6 +55,7 @@ module mld_z_onelev_mod use mld_base_prec_type use mld_z_base_smoother_mod + use mld_z_base_aggregator_mod use psb_base_mod, only : psb_zspmat_type, psb_z_vect_type, & & psb_z_base_vect_type, psb_zlinmap_type, psb_dpk_, & & psb_ipk_, psb_long_int_k_, psb_desc_type, psb_i_base_vect_type, & @@ -136,9 +137,10 @@ module mld_z_onelev_mod & z_wrk_clone, z_wrk_move_alloc, z_wrk_cnv type mld_z_onelev_type - class(mld_z_base_smoother_type), allocatable :: sm, sm2a - class(mld_z_base_smoother_type), pointer :: sm2 => null() - class(mld_zmlprec_wrk_type), allocatable :: wrk + class(mld_z_base_smoother_type), allocatable :: sm, sm2a + class(mld_z_base_smoother_type), pointer :: sm2 => null() + class(mld_zmlprec_wrk_type), allocatable :: wrk + class(mld_z_base_aggregator_type), allocatable :: aggr type(mld_dml_parms) :: parms type(psb_zspmat_type) :: ac integer(psb_ipk_) :: ac_nz_loc, ac_nz_tot @@ -166,8 +168,9 @@ module mld_z_onelev_mod procedure, pass(lv) :: csetc => mld_z_base_onelev_csetc procedure, pass(lv) :: setsm => mld_z_base_onelev_setsm procedure, pass(lv) :: setsv => mld_z_base_onelev_setsv + procedure, pass(lv) :: setag => mld_z_base_onelev_setag generic, public :: set => seti, setr, setc, & - & cseti, csetr, csetc, setsm, setsv + & cseti, csetr, csetc, setsm, setsv, setag procedure, pass(lv) :: sizeof => z_base_onelev_sizeof procedure, pass(lv) :: get_nzeros => z_base_onelev_get_nzeros procedure, pass(lv) :: get_wrksz => z_base_onelev_get_wrksize @@ -299,6 +302,20 @@ module mld_z_onelev_mod end subroutine mld_z_base_onelev_setsv end interface + interface + subroutine mld_z_base_onelev_setag(lv,val,info,pos) + import :: psb_dpk_, mld_z_onelev_type, mld_z_base_aggregator_type, & + & psb_ipk_, psb_long_int_k_, psb_desc_type + Implicit None + + ! Arguments + class(mld_z_onelev_type), target, intent(inout) :: lv + class(mld_z_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + end subroutine mld_z_base_onelev_setag + end interface + interface subroutine mld_z_base_onelev_setc(lv,what,val,info,pos) import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & @@ -323,9 +340,9 @@ module mld_z_onelev_mod class(mld_z_onelev_type), intent(inout) :: lv integer(psb_ipk_), intent(in) :: what - real(psb_dpk_), intent(in) :: val + real(psb_dpk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info - character(len=*), optional, intent(in) :: pos + character(len=*), optional, intent(in) :: pos end subroutine mld_z_base_onelev_setr end interface @@ -450,7 +467,8 @@ contains Implicit None ! Arguments - class(mld_z_onelev_type), target, intent(inout) :: lv + class(mld_z_onelev_type), target, intent(inout) :: lv + integer(psb_ipk_) :: info lv%parms%sweeps_pre = 1 lv%parms%sweeps_post = 1 @@ -473,7 +491,9 @@ contains else lv%sm2 => lv%sm end if - + if (.not.allocated(lv%aggr)) allocate(mld_z_base_aggregator_type :: lv%aggr,stat=info) + if (allocated(lv%aggr)) call lv%aggr%default() + return end subroutine z_base_onelev_default @@ -508,6 +528,14 @@ contains end if lvout%sm2 => lvout%sm end if + if (allocated(lv%aggr)) then + call lv%aggr%clone(lvout%aggr,info) + else + if (allocated(lvout%aggr)) then + call lvout%aggr%free(info) + if (info==psb_success_) deallocate(lvout%aggr,stat=info) + end if + end if if (info == psb_success_) call lv%parms%clone(lvout%parms,info) if (info == psb_success_) call lv%ac%clone(lvout%ac,info) if (info == psb_success_) call lv%tprol%clone(lvout%tprol,info) @@ -538,8 +566,9 @@ contains call move_alloc(lv%sm2a,b%sm2a) b%sm2 =>b%sm end if - - if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) + + call move_alloc(lv%aggr,b%aggr) + if (info == psb_success_) call psb_move_alloc(lv%ac,b%ac,info) if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(lv%map,b%map,info) diff --git a/mlprec/mld_z_prec_mod.f90 b/mlprec/mld_z_prec_mod.f90 index 1c094a84..374a41b2 100644 --- a/mlprec/mld_z_prec_mod.f90 +++ b/mlprec/mld_z_prec_mod.f90 @@ -55,7 +55,8 @@ module mld_z_prec_mod interface mld_precset module procedure mld_z_iprecsetsm, mld_z_iprecsetsv, & & mld_z_iprecseti, mld_z_iprecsetc, mld_z_iprecsetr, & - & mld_z_cprecseti, mld_z_cprecsetc, mld_z_cprecsetr + & mld_z_cprecseti, mld_z_cprecsetc, mld_z_cprecsetr, & + & mld_z_iprecsetag end interface mld_precset interface mld_extprol_bld @@ -97,6 +98,14 @@ contains call p%set(val,info, pos=pos) end subroutine mld_z_iprecsetsv + subroutine mld_z_iprecsetag(p,val,info,pos) + type(mld_zprec_type), intent(inout) :: p + class(mld_z_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + character(len=*), optional, intent(in) :: pos + call p%set(val,info, pos=pos) + end subroutine mld_z_iprecsetag + subroutine mld_z_iprecseti(p,what,val,info,pos) type(mld_zprec_type), intent(inout) :: p integer(psb_ipk_), intent(in) :: what diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 index 39b11f96..3def4f30 100644 --- a/mlprec/mld_z_prec_type.f90 +++ b/mlprec/mld_z_prec_type.f90 @@ -55,6 +55,7 @@ module mld_z_prec_type use mld_base_prec_type use mld_z_base_solver_mod use mld_z_base_smoother_mod + use mld_z_base_aggregator_mod use mld_z_onelev_mod use psb_prec_mod, only : psb_zprec_type @@ -92,8 +93,8 @@ module mld_z_prec_type ! 2. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_levs = 20_psb_ipk_ ! 3. min_cr_ratio = 1.5 - real(psb_dpk_) :: min_cr_ratio = 1.5_psb_dpk_ - real(psb_dpk_) :: op_complexity=dzero + real(psb_dpk_) :: min_cr_ratio = 1.5_psb_dpk_ + real(psb_dpk_) :: op_complexity = dzero ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -126,6 +127,7 @@ module mld_z_prec_type procedure, pass(prec) :: sizeof => mld_zprec_sizeof procedure, pass(prec) :: setsm => mld_zprecsetsm procedure, pass(prec) :: setsv => mld_zprecsetsv + procedure, pass(prec) :: setag => mld_zprecsetag procedure, pass(prec) :: seti => mld_zprecseti procedure, pass(prec) :: setc => mld_zprecsetc procedure, pass(prec) :: setr => mld_zprecsetr @@ -133,7 +135,7 @@ module mld_z_prec_type procedure, pass(prec) :: csetc => mld_zcprecsetc procedure, pass(prec) :: csetr => mld_zcprecsetr generic, public :: set => seti, setc, setr, & - & cseti, csetc, csetr, setsm, setsv + & cseti, csetc, csetr, setsm, setsv, setag procedure, pass(prec) :: get_smoother => mld_z_get_smootherp procedure, pass(prec) :: get_solver => mld_z_get_solverp procedure, pass(prec) :: move_alloc => z_prec_move_alloc @@ -234,6 +236,15 @@ module mld_z_prec_type integer(psb_ipk_), optional, intent(in) :: ilev,ilmax character(len=*), optional, intent(in) :: pos end subroutine mld_zprecsetsv + subroutine mld_zprecsetag(prec,val,info,ilev,pos) + import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & + & mld_zprec_type, mld_z_base_aggregator_type, psb_ipk_ + class(mld_zprec_type), intent(inout) :: prec + class(mld_z_base_aggregator_type), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: ilev + character(len=*), optional, intent(in) :: pos + end subroutine mld_zprecsetag subroutine mld_zprecseti(prec,what,val,info,ilev,ilmax,pos) import :: psb_zspmat_type, psb_desc_type, psb_dpk_, & & mld_zprec_type, psb_ipk_ diff --git a/mlprec/mld_z_symdec_aggregator_mod.f90 b/mlprec/mld_z_symdec_aggregator_mod.f90 new file mode 100644 index 00000000..f1f08e7d --- /dev/null +++ b/mlprec/mld_z_symdec_aggregator_mod.f90 @@ -0,0 +1,128 @@ +! +! +! MLD2P4 version 2.1 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008, 2010, 2012, 2015, 2017 , 2017 +! +! Salvatore Filippone Cranfield University +! Ambra Abdullahi Hassan University of Rome Tor Vergata +! Pasqua D'Ambra IAC-CNR, Naples, IT +! Daniela di Serafino University of Campania "L. Vanvitelli", Caserta, IT +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the MLD2P4 group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This version only differs from the +! basic decoupled aggregation algorithm because it works on (the +! pattern of) A+A^T instead of A. +! +! M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! P. D'Ambra, D. di Serafino and S. Filippone, On the development of +! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. +! 57 (2007), 1181-1196. +! +module mld_z_symdec_aggregator_mod + + use mld_z_base_aggregator_mod + ! + ! sm - class(mld_T_base_smoother_type), allocatable + ! The current level preconditioner (aka smoother). + ! parms - type(mld_RTml_parms) + ! The parameters defining the multilevel strategy. + ! ac - The local part of the current-level matrix, built by + ! coarsening the previous-level matrix. + ! desc_ac - type(psb_desc_type). + ! The communication descriptor associated to the matrix + ! stored in ac. + ! base_a - type(psb_Tspmat_type), pointer. + ! Pointer (really a pointer!) to the local part of the current + ! matrix (so we have a unified treatment of residuals). + ! We need this to avoid passing explicitly the current matrix + ! to the routine which applies the preconditioner. + ! base_desc - type(psb_desc_type), pointer. + ! Pointer to the communication descriptor associated to the + ! matrix pointed by base_a. + ! map - Stores the maps (restriction and prolongation) between the + ! vector spaces associated to the index spaces of the previous + ! and current levels. + ! + ! Methods: + ! Most methods follow the encapsulation hierarchy: they take whatever action + ! is appropriate for the current object, then call the corresponding method for + ! the contained object. + ! As an example: the descr() method prints out a description of the + ! level. It starts by invoking the descr() method of the parms object, + ! then calls the descr() method of the smoother object. + ! + ! descr - Prints a description of the object. + ! default - Set default values + ! dump - Dump to file object contents + ! set - Sets various parameters; when a request is unknown + ! it is passed to the smoother object for further processing. + ! check - Sanity checks. + ! sizeof - Total memory occupation in bytes + ! get_nzeros - Number of nonzeros + ! + ! + type, extends(mld_z_base_aggregator_type) :: mld_z_symdec_aggregator_type + + contains + procedure, pass(ag) :: tprol => mld_z_symdec_aggregator_build_tprol + procedure, nopass :: fmt => mld_z_symdec_aggregator_fmt + end type mld_z_symdec_aggregator_type + + + interface + subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + import :: mld_z_symdec_aggregator_type, psb_desc_type, psb_zspmat_type, psb_dpk_, & + & psb_ipk_, psb_long_int_k_, mld_dml_parms + implicit none + class(mld_z_symdec_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_symdec_aggregator_build_tprol + end interface + + +contains + + function mld_z_symdec_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Symmetric Decoupled aggregation" + end function mld_z_symdec_aggregator_fmt + +end module mld_z_symdec_aggregator_mod