diff --git a/mlprec/impl/Makefile b/mlprec/impl/Makefile index a2c35c31..9a0d92bf 100644 --- a/mlprec/impl/Makefile +++ b/mlprec/impl/Makefile @@ -8,13 +8,13 @@ FINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(MODDIR) $(FMFLAG)$(INCDIR) $(PSBLAS_INCLUD -DMPFOBJS=mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o mld_daggrmat_biz_asb.o +DMPFOBJS= -SMPFOBJS=mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o mld_saggrmat_minnrg_asb.o mld_saggrmat_biz_asb.o +SMPFOBJS= -ZMPFOBJS=mld_zaggrmat_nosmth_asb.o mld_zaggrmat_smth_asb.o mld_zaggrmat_minnrg_asb.o mld_zaggrmat_biz_asb.o +ZMPFOBJS= -CMPFOBJS=mld_caggrmat_nosmth_asb.o mld_caggrmat_smth_asb.o mld_caggrmat_minnrg_asb.o mld_caggrmat_biz_asb.o +CMPFOBJS= MPFOBJS=$(SMPFOBJS) $(DMPFOBJS) $(CMPFOBJS) $(ZMPFOBJS) @@ -24,26 +24,22 @@ MPCOBJS=mld_dslud_interface.o mld_zslud_interface.o DINNEROBJS= mld_dmlprec_bld.o mld_dfile_prec_descr.o \ mld_d_smoothers_bld.o mld_d_hierarchy_bld.o \ - mld_daggrmap_bld.o \ - mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \ - $(DMPFOBJS) mld_d_extprol_bld.o mld_d_lev_aggrmap_bld.o mld_d_lev_aggrmat_asb.o + mld_dmlprec_aply.o \ + $(DMPFOBJS) mld_d_extprol_bld.o mld_d_lev_aggrmap_bld.o mld_d_lev_aggrmat_asb.o SINNEROBJS= mld_smlprec_bld.o mld_sfile_prec_descr.o \ mld_s_smoothers_bld.o mld_s_hierarchy_bld.o \ - mld_saggrmap_bld.o \ - mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \ + mld_smlprec_aply.o \ $(SMPFOBJS) mld_s_extprol_bld.o mld_s_lev_aggrmap_bld.o mld_s_lev_aggrmat_asb.o ZINNEROBJS= mld_zmlprec_bld.o mld_zfile_prec_descr.o \ mld_z_smoothers_bld.o mld_z_hierarchy_bld.o \ - mld_zaggrmap_bld.o \ - mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \ + mld_zmlprec_aply.o \ $(ZMPFOBJS) mld_z_extprol_bld.o mld_z_lev_aggrmap_bld.o mld_z_lev_aggrmat_asb.o CINNEROBJS= mld_cmlprec_bld.o mld_cfile_prec_descr.o \ mld_c_smoothers_bld.o mld_c_hierarchy_bld.o \ - mld_caggrmap_bld.o \ - mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \ + mld_cmlprec_aply.o \ $(CMPFOBJS) mld_c_extprol_bld.o mld_c_lev_aggrmap_bld.o mld_c_lev_aggrmat_asb.o INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS) @@ -71,10 +67,12 @@ OBJS=$(F90OBJS) $(COBJS) $(MPCOBJS) LIBNAME=libmld_prec.a -lib: $(OBJS) levd smoothd solvd +lib: $(OBJS) levd smoothd solvd aggrd $(AR) $(HERE)/$(LIBNAME) $(OBJS) $(RANLIB) $(HERE)/$(LIBNAME) +aggrd: + cd aggregator && $(MAKE) levd: cd level && $(MAKE) smoothd: @@ -89,9 +87,11 @@ mpobjs: veryclean: clean /bin/rm -f $(LIBNAME) -clean: solvclean smoothclean levclean +clean: solvclean smoothclean levclean aggrclean /bin/rm -f $(OBJS) $(LOCAL_MODS) +aggrclean: + cd aggregator && $(MAKE) clean levclean: cd level && $(MAKE) clean smoothclean: diff --git a/mlprec/impl/aggregator/Makefile b/mlprec/impl/aggregator/Makefile new file mode 100644 index 00000000..762e09eb --- /dev/null +++ b/mlprec/impl/aggregator/Makefile @@ -0,0 +1,61 @@ +include ../../../Make.inc +LIBDIR=../../../lib +INCDIR=../../../include +MODDIR=../../../modules +HERE=../.. + +FINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(MODDIR) $(FMFLAG)$(INCDIR) $(PSBLAS_INCLUDES) + +#CINCLUDES= -I${SUPERLU_INCDIR} -I${HSL_INCDIR} -I${SPRAL_INCDIR} -I/home/users/pasqua/Ambra/BootCMatch/include -lBCM -L/home/users/pasqua/Ambra/BootCMatch/lib -lm + +OBJS= \ +bootCMatch_interface.o\ +mld_s_base_aggregator_mat_asb.o \ +mld_s_base_aggregator_tprol.o \ +mld_s_hybrid_aggregator_tprol.o \ +mld_s_symdec_aggregator_tprol.o \ +mld_s_map_to_tprol.o mld_s_dec_map_bld.o mld_s_hyb_map_bld.o\ +mld_saggrmat_biz_asb.o mld_saggrmat_minnrg_asb.o\ +mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o \ +mld_d_base_aggregator_mat_asb.o \ +mld_d_base_aggregator_tprol.o \ +mld_d_bcmatch_aggregator_tprol.o\ +mld_d_hybrid_aggregator_tprol.o \ +mld_d_symdec_aggregator_tprol.o \ +mld_d_map_to_tprol.o mld_d_dec_map_bld.o mld_d_hyb_map_bld.o\ +mld_daggrmat_unsmth_spmm_asb.o\ +mld_daggrmat_biz_asb.o mld_daggrmat_minnrg_asb.o\ +mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o \ +mld_c_base_aggregator_mat_asb.o \ +mld_c_base_aggregator_tprol.o \ +mld_c_hybrid_aggregator_tprol.o \ +mld_c_symdec_aggregator_tprol.o \ +mld_c_map_to_tprol.o mld_c_dec_map_bld.o mld_c_hyb_map_bld.o\ +mld_caggrmat_biz_asb.o mld_caggrmat_minnrg_asb.o\ +mld_caggrmat_nosmth_asb.o mld_caggrmat_smth_asb.o \ +mld_z_base_aggregator_mat_asb.o \ +mld_z_base_aggregator_tprol.o \ +mld_z_hybrid_aggregator_tprol.o \ +mld_z_symdec_aggregator_tprol.o \ +mld_z_map_to_tprol.o mld_z_dec_map_bld.o mld_z_hyb_map_bld.o\ +mld_zaggrmat_biz_asb.o mld_zaggrmat_minnrg_asb.o\ +mld_zaggrmat_nosmth_asb.o mld_zaggrmat_smth_asb.o + +#mld_d_bcmatch_map_to_tprol.o mld_d_bcmatch_aggregator_mat_asb.o \ + + +LIBNAME=libmld_prec.a + +lib: $(OBJS) + $(AR) $(HERE)/$(LIBNAME) $(OBJS) + $(RANLIB) $(HERE)/$(LIBNAME) + +mpobjs: + (make $(MPFOBJS) F90="$(MPF90)" F90COPT="$(F90COPT)") + (make $(MPCOBJS) CC="$(MPCC)" CCOPT="$(CCOPT)") + +veryclean: clean + /bin/rm -f $(LIBNAME) + +clean: + /bin/rm -f $(OBJS) $(LOCAL_MODS) diff --git a/mlprec/impl/aggregator/mld_c_base_aggregator_mat_asb.f90 b/mlprec/impl/aggregator/mld_c_base_aggregator_mat_asb.f90 new file mode 100644 index 00000000..dc4bb1df --- /dev/null +++ b/mlprec/impl/aggregator/mld_c_base_aggregator_mat_asb.f90 @@ -0,0 +1,225 @@ +! +! +! 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. +! +! +! File: mld_c_base_aggregator_mat_asb.f90 +! +! Subroutine: mld_c_base_aggregator_mat_asb +! Version: complex +! +! This routine builds the matrix associated to the current level of the +! multilevel preconditioner from the matrix associated to the previous level, +! by using the user-specified aggregation technique (therefore, it also builds the +! prolongation and restriction operators mapping the current level to the +! previous one and vice versa). +! The current level is regarded as the coarse one, while the previous as +! the fine one. This is in agreement with the fact that the routine is called, +! by mld_mlprec_bld, only on levels >=2. +! The coarse-level matrix A_C is built from a fine-level matrix A +! by using the Galerkin approach, i.e. +! +! A_C = P_C^T A P_C, +! +! where P_C is a prolongator from the coarse level to the fine one. +! +! A mapping from the nodes of the adjacency graph of A to the nodes of the +! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine. +! The prolongator P_C is built here from this mapping, according to the +! value of p%iprcparm(mld_aggr_kind_), specified by the user through +! mld_cprecinit and mld_zprecset. +! On output from this routine the entries of AC, op_prol, op_restr +! are still in "global numbering" mode; this is fixed in the calling routine +! mld_c_lev_aggrmat_asb. +! +! Currently four different prolongators are implemented, corresponding to +! four aggregation algorithms: +! 1. un-smoothed aggregation, +! 2. smoothed aggregation, +! 3. "bizarre" aggregation. +! 4. minimum energy +! 1. The non-smoothed aggregation uses as prolongator the piecewise constant +! interpolation operator corresponding to the fine-to-coarse level mapping built +! by p%aggr%bld_tprol. This is called tentative prolongator. +! 2. The smoothed aggregation uses as prolongator the operator obtained by applying +! a damped Jacobi smoother to the tentative prolongator. +! 3. The "bizarre" aggregation uses a prolongator proposed by the authors of MLD2P4. +! This prolongator still requires a deep analysis and testing and its use is +! not recommended. +! 4. Minimum energy aggregation +! +! For more details see +! 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. +! M. Sala, R. Tuminaro: A new Petrov-Galerkin smoothed aggregation preconditioner +! for nonsymmetric linear systems, SIAM J. Sci. Comput., 31(1):143-166 (2008) +! +! +! The main structure is: +! 1. Perform sanity checks; +! 2. Compute prolongator/restrictor/AC +! +! +! Arguments: +! ag - type(mld_c_base_aggregator_type), input/output. +! The aggregator object +! parms - type(mld_sml_parms), input +! The aggregation parameters +! a - type(psb_cspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! The 'one-level' data structure that will contain the local +! part of the matrix to be built as well as the information +! concerning the prolongator and its transpose. +! ilaggr - integer, dimension(:), input +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that the indices +! are assumed to be shifted so as to make sure the ranges on +! the various processes do not overlap. +! nlaggr - integer, dimension(:) input +! nlaggr(i) contains the aggregates held by process i. +! ac - type(psb_cspmat_type), output +! The coarse matrix on output +! +! op_prol - type(psb_cspmat_type), input/output +! The tentative prolongator on input, the computed prolongator on output +! +! op_restr - type(psb_cspmat_type), output +! The restrictor operator; normally, it is the transpose of the prolongator. +! +! info - integer, output. +! Error code. +! +subroutine mld_c_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info) + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_base_aggregator_mat_asb + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + type(psb_c_coo_sparse_mat) :: acoo, bcoo + type(psb_c_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl,ntaggr + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_c_base_aggregator_mat_asb' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%aggr_kind,'Smoother',& + & mld_smooth_prol_,is_legal_ml_aggr_kind) + call mld_check_def(parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + call mld_check_def(parms%aggr_filter,'Use filtered matrix',& + & mld_no_filter_mat_,is_legal_aggr_filter) + call mld_check_def(parms%smoother_pos,'smooth_pos',& + & mld_pre_smooth_,is_legal_ml_smooth_pos) + call mld_check_def(parms%aggr_omega_alg,'Omega Alg.',& + & mld_eig_est_,is_legal_ml_aggr_omega_alg) + call mld_check_def(parms%aggr_eig,'Eigenvalue estimate',& + & mld_max_norm_,is_legal_ml_aggr_eig) + call mld_check_def(parms%aggr_omega_val,'Omega',szero,is_legal_s_omega) + + ! + ! Build the coarse-level matrix from the fine-level one, starting from + ! the mapping defined by mld_aggrmap_bld and applying the aggregation + ! algorithm specified by p%iprcparm(mld_aggr_kind_) + ! + select case (parms%aggr_kind) + case (mld_no_smooth_) + + call mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & parms,ac,op_prol,op_restr,info) + + case(mld_smooth_prol_) + + call mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_biz_prol_) + + call mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + +end subroutine mld_c_base_aggregator_mat_asb diff --git a/mlprec/impl/aggregator/mld_c_base_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_c_base_aggregator_tprol.f90 new file mode 100644 index 00000000..18c5428d --- /dev/null +++ b/mlprec/impl/aggregator/mld_c_base_aggregator_tprol.f90 @@ -0,0 +1,124 @@ +! +! +! 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. +! +! +! File: mld_c_base_aggregator_tprol.f90 +! +! Subroutine: mld_c_base_aggregator_tprol +! Version: complex +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_c_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_cspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_cspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_c_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod +! use mld_c_base_aggregator_mod + use mld_c_inner_mod, mld_protect_name => mld_c_base_aggregator_build_tprol + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_c_base_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) + + + call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_c_base_aggregator_build_tprol diff --git a/mlprec/impl/mld_c_dec_map_bld.f90 b/mlprec/impl/aggregator/mld_c_dec_map_bld.f90 similarity index 100% rename from mlprec/impl/mld_c_dec_map_bld.f90 rename to mlprec/impl/aggregator/mld_c_dec_map_bld.f90 diff --git a/mlprec/impl/aggregator/mld_c_hyb_map_bld.f90 b/mlprec/impl/aggregator/mld_c_hyb_map_bld.f90 new file mode 100644 index 00000000..22ad0a2c --- /dev/null +++ b/mlprec/impl/aggregator/mld_c_hyb_map_bld.f90 @@ -0,0 +1,328 @@ +! +! +! 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. +! +! +! +! File: mld_c_hyb_map__bld.f90 +! +! Subroutine: mld_c_hyb_map_bld +! Version: complex +! +! This routine builds the tentative prolongator based on 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. +! +! Note: upon exit +! +! Arguments: +! a - type(psb_cspmat_type). +! The sparse matrix structure containing the local part of the +! matrix to be preconditioned. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_cprec_type), input/output. +! The preconditioner data structure; upon exit it contains +! the multilevel hierarchy of prolongators, restrictors +! and coarse matrices. +! info - integer, output. +! Error code. +! +! +! +subroutine mld_c_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_hyb_map_bld + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:), tmpaggr(:) + complex(psb_spk_), allocatable :: val(:), diag(:) + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt + type(psb_c_csr_sparse_mat) :: acsr, muij, s_neigh + type(psb_c_coo_sparse_mat) :: s_neigh_coo + real(psb_spk_) :: cpling, tcl + logical :: disjoint + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if (psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_hyb_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + nr = a%get_nrows() + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),icol(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + diag = a%get_diag(info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + ! + ! Phase zero: compute muij + ! + call a%cp_to(muij) + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) + end do + end do + !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! + ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 + ! + call s_neigh_coo%allocate(nr,nr,muij%get_nzeros()) + ip = 0 + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + ip = ip + 1 + s_neigh_coo%ia(ip) = i + s_neigh_coo%ja(ip) = j + if (real(muij%val(k)) >= theta) then + s_neigh_coo%val(ip) = sone + else + s_neigh_coo%val(ip) = -sone + end if + end do + end do + !write(*,*) 'S_NEIGH: ',nr,ip + call s_neigh_coo%set_nzeros(ip) + call s_neigh%mv_from_coo(s_neigh_coo,info) + + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = muij%irp(i+1) - muij%irp(i) + end do + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if + + + ! + ! Phase one: Start with disjoint groups. + ! + naggr = 0 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Get the 1-neighbourhood of I + ! + ip1 = s_neigh%irp(i) + nz = s_neigh%irp(i+1)-ip1 + ! + ! If the neighbourhood only contains I, skip it + ! + if (nz ==0) then + ilaggr(i) = 0 + cycle step1 + end if + if ((nz==1).and.(s_neigh%ja(ip1)==i)) then + ilaggr(i) = 0 + cycle step1 + end if + ! + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate. + ! + nzcnt = count(real(s_neigh%val(ip1:ip1+nz-1)) > 0) + icol(1:nzcnt) = pack(s_neigh%ja(ip1:ip1+nz-1),(real(s_neigh%val(ip1:ip1+nz-1)) > 0)) + disjoint = all(ilaggr(icol(1:nzcnt)) == -(nr+1)) + if (disjoint) then + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, nzcnt + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 + + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: join the neighbours + ! + tmpaggr = ilaggr + step2: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = szero + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if ( (tmpaggr(j) > 0).and. (real(muij%val(k)) > cpling)& + & .and.(real(s_neigh%val(k))>0)) then + ip = k + cpling = muij%val(k) + end if + end if + enddo + if (ip > 0) then + ilaggr(i) = ilaggr(s_neigh%ja(ip)) + end if + end if + end do step2 + + + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) < 0) then + ! + ! Find its strongly connected neighbourhood not + ! already aggregated, and make it into a new aggregate. + ! + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if (ilaggr(j) < 0) then + ip = ip + 1 + icol(ip) = j + end if + end if + enddo + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do + end if + end if + end do step3 + + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_hyb_map_bld + diff --git a/mlprec/impl/aggregator/mld_c_hybrid_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_c_hybrid_aggregator_tprol.f90 new file mode 100644 index 00000000..675baaa7 --- /dev/null +++ b/mlprec/impl/aggregator/mld_c_hybrid_aggregator_tprol.f90 @@ -0,0 +1,122 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_c_hybrid_aggregator_tprol.f90 +! +! Subroutine: mld_c_hybrid_aggregator_tprol +! Version: complex +! +! +! This routine is mainly an interface to hyb_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_c_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_cspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_cspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_c_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_c_hybrid_aggregator_mod, mld_protect_name => mld_c_hybrid_aggregator_build_tprol + use mld_c_inner_mod + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_c_hybrid_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(p%parms%ml_cycle,'Multilevel cycle',& + & mld_mult_ml_,is_legal_ml_cycle) + call mld_check_def(p%parms%par_aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_par_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) + call mld_hyb_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_c_hybrid_aggregator_build_tprol diff --git a/mlprec/impl/aggregator/mld_c_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_c_map_to_tprol.f90 new file mode 100644 index 00000000..ae2154c5 --- /dev/null +++ b/mlprec/impl/aggregator/mld_c_map_to_tprol.f90 @@ -0,0 +1,150 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_c_map_to_tprol.f90 +! +! Subroutine: mld_c_map_to_tprol +! Version: complex +! +! This routine uses a mapping from the row indices of the fine-level matrix +! to the row indices of the coarse-level matrix to build a tentative +! prolongator, i.e. a piecewise constant operator. +! This is later used to build the final operator; the code has been refactored here +! to be shared among all the methods that provide the tentative prolongator +! through a simple integer mapping. +! +! The aggregation algorithm is a parallel version of that described in +! * M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! * P. Vanek, J. Mandel and M. Brezina, Algebraic Multigrid by Smoothed +! Aggregation for Second and Fourth Order Elliptic Problems, Computing, 56 +! (1996), 179-196. +! For more details see +! 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. +! +! +! Arguments: +! aggr_type - integer, input. +! The scalar used to identify the aggregation algorithm. +! theta - real, input. +! The aggregation threshold used in the aggregation algorithm. +! a - type(psb_cspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! ilaggr - integer, dimension(:), allocatable. +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable. +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_cspmat_type). +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_c_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_map_to_tprol + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr + type(psb_c_coo_sparse_mat) :: tmpcoo + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_map_to_tprol' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 + call psb_halo(ilaggr,desc_a,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') + goto 9999 + end if + + call tmpcoo%allocate(ncol,ntaggr,ncol) + do i=1,ncol + tmpcoo%val(i) = cone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) + end do + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_sorted() ! At this point this is in row-major + call op_prol%mv_from(tmpcoo) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_map_to_tprol diff --git a/mlprec/impl/aggregator/mld_c_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_c_symdec_aggregator_tprol.f90 new file mode 100644 index 00000000..8fe0b7a8 --- /dev/null +++ b/mlprec/impl/aggregator/mld_c_symdec_aggregator_tprol.f90 @@ -0,0 +1,141 @@ +! +! +! 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. +! +! +! File: mld_c_symdec_aggregator_tprol.f90 +! +! Subroutine: mld_c_symdec_aggregator_tprol +! Version: complex +! +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! +! Arguments: +! p - type(mld_c_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_cspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_cspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_c_symdec_aggregator_mod, mld_protect_name => mld_c_symdec_aggregator_build_tprol + use mld_c_inner_mod + 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 + + ! Local variables + type(psb_cspmat_type) :: atmp, atrans + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr, nr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_c_symdec_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) + + nr = a%get_nrows() + call a%csclip(atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atmp%transp(atrans) + if (info == psb_success_) call atrans%cscnv(info,type='COO') + if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atrans%free() + if (info == psb_success_) call atmp%cscnv(info,type='CSR') + + if (info == psb_success_) & + & call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call atmp%free() + + if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_c_symdec_aggregator_build_tprol diff --git a/mlprec/impl/mld_caggrmat_biz_asb.f90 b/mlprec/impl/aggregator/mld_caggrmat_biz_asb.f90 similarity index 100% rename from mlprec/impl/mld_caggrmat_biz_asb.f90 rename to mlprec/impl/aggregator/mld_caggrmat_biz_asb.f90 diff --git a/mlprec/impl/mld_caggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_caggrmat_minnrg_asb.f90 similarity index 100% rename from mlprec/impl/mld_caggrmat_minnrg_asb.f90 rename to mlprec/impl/aggregator/mld_caggrmat_minnrg_asb.f90 diff --git a/mlprec/impl/mld_caggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_caggrmat_nosmth_asb.f90 similarity index 100% rename from mlprec/impl/mld_caggrmat_nosmth_asb.f90 rename to mlprec/impl/aggregator/mld_caggrmat_nosmth_asb.f90 diff --git a/mlprec/impl/mld_caggrmat_smth_asb.f90 b/mlprec/impl/aggregator/mld_caggrmat_smth_asb.f90 similarity index 100% rename from mlprec/impl/mld_caggrmat_smth_asb.f90 rename to mlprec/impl/aggregator/mld_caggrmat_smth_asb.f90 diff --git a/mlprec/impl/aggregator/mld_d_base_aggregator_mat_asb.f90 b/mlprec/impl/aggregator/mld_d_base_aggregator_mat_asb.f90 new file mode 100644 index 00000000..ba600776 --- /dev/null +++ b/mlprec/impl/aggregator/mld_d_base_aggregator_mat_asb.f90 @@ -0,0 +1,225 @@ +! +! +! 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. +! +! +! File: mld_d_base_aggregator_mat_asb.f90 +! +! Subroutine: mld_d_base_aggregator_mat_asb +! Version: real +! +! This routine builds the matrix associated to the current level of the +! multilevel preconditioner from the matrix associated to the previous level, +! by using the user-specified aggregation technique (therefore, it also builds the +! prolongation and restriction operators mapping the current level to the +! previous one and vice versa). +! The current level is regarded as the coarse one, while the previous as +! the fine one. This is in agreement with the fact that the routine is called, +! by mld_mlprec_bld, only on levels >=2. +! The coarse-level matrix A_C is built from a fine-level matrix A +! by using the Galerkin approach, i.e. +! +! A_C = P_C^T A P_C, +! +! where P_C is a prolongator from the coarse level to the fine one. +! +! A mapping from the nodes of the adjacency graph of A to the nodes of the +! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine. +! The prolongator P_C is built here from this mapping, according to the +! value of p%iprcparm(mld_aggr_kind_), specified by the user through +! mld_dprecinit and mld_zprecset. +! On output from this routine the entries of AC, op_prol, op_restr +! are still in "global numbering" mode; this is fixed in the calling routine +! mld_d_lev_aggrmat_asb. +! +! Currently four different prolongators are implemented, corresponding to +! four aggregation algorithms: +! 1. un-smoothed aggregation, +! 2. smoothed aggregation, +! 3. "bizarre" aggregation. +! 4. minimum energy +! 1. The non-smoothed aggregation uses as prolongator the piecewise constant +! interpolation operator corresponding to the fine-to-coarse level mapping built +! by p%aggr%bld_tprol. This is called tentative prolongator. +! 2. The smoothed aggregation uses as prolongator the operator obtained by applying +! a damped Jacobi smoother to the tentative prolongator. +! 3. The "bizarre" aggregation uses a prolongator proposed by the authors of MLD2P4. +! This prolongator still requires a deep analysis and testing and its use is +! not recommended. +! 4. Minimum energy aggregation +! +! For more details see +! 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. +! M. Sala, R. Tuminaro: A new Petrov-Galerkin smoothed aggregation preconditioner +! for nonsymmetric linear systems, SIAM J. Sci. Comput., 31(1):143-166 (2008) +! +! +! The main structure is: +! 1. Perform sanity checks; +! 2. Compute prolongator/restrictor/AC +! +! +! Arguments: +! ag - type(mld_d_base_aggregator_type), input/output. +! The aggregator object +! parms - type(mld_dml_parms), input +! The aggregation parameters +! a - type(psb_dspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! The 'one-level' data structure that will contain the local +! part of the matrix to be built as well as the information +! concerning the prolongator and its transpose. +! ilaggr - integer, dimension(:), input +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that the indices +! are assumed to be shifted so as to make sure the ranges on +! the various processes do not overlap. +! nlaggr - integer, dimension(:) input +! nlaggr(i) contains the aggregates held by process i. +! ac - type(psb_dspmat_type), output +! The coarse matrix on output +! +! op_prol - type(psb_dspmat_type), input/output +! The tentative prolongator on input, the computed prolongator on output +! +! op_restr - type(psb_dspmat_type), output +! The restrictor operator; normally, it is the transpose of the prolongator. +! +! info - integer, output. +! Error code. +! +subroutine mld_d_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info) + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_base_aggregator_mat_asb + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + type(psb_d_coo_sparse_mat) :: acoo, bcoo + type(psb_d_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl,ntaggr + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_d_base_aggregator_mat_asb' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%aggr_kind,'Smoother',& + & mld_smooth_prol_,is_legal_ml_aggr_kind) + call mld_check_def(parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + call mld_check_def(parms%aggr_filter,'Use filtered matrix',& + & mld_no_filter_mat_,is_legal_aggr_filter) + call mld_check_def(parms%smoother_pos,'smooth_pos',& + & mld_pre_smooth_,is_legal_ml_smooth_pos) + call mld_check_def(parms%aggr_omega_alg,'Omega Alg.',& + & mld_eig_est_,is_legal_ml_aggr_omega_alg) + call mld_check_def(parms%aggr_eig,'Eigenvalue estimate',& + & mld_max_norm_,is_legal_ml_aggr_eig) + call mld_check_def(parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega) + + ! + ! Build the coarse-level matrix from the fine-level one, starting from + ! the mapping defined by mld_aggrmap_bld and applying the aggregation + ! algorithm specified by p%iprcparm(mld_aggr_kind_) + ! + select case (parms%aggr_kind) + case (mld_no_smooth_) + + call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & parms,ac,op_prol,op_restr,info) + + case(mld_smooth_prol_) + + call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_biz_prol_) + + call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + +end subroutine mld_d_base_aggregator_mat_asb diff --git a/mlprec/impl/aggregator/mld_d_base_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_d_base_aggregator_tprol.f90 new file mode 100644 index 00000000..02ffb7d7 --- /dev/null +++ b/mlprec/impl/aggregator/mld_d_base_aggregator_tprol.f90 @@ -0,0 +1,124 @@ +! +! +! 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. +! +! +! File: mld_d_base_aggregator_tprol.f90 +! +! Subroutine: mld_d_base_aggregator_tprol +! Version: real +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_d_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_dspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_dspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_d_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod +! use mld_d_base_aggregator_mod + use mld_d_inner_mod, mld_protect_name => mld_d_base_aggregator_build_tprol + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_d_base_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + + + call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_d_base_aggregator_build_tprol diff --git a/mlprec/impl/mld_d_dec_map_bld.f90 b/mlprec/impl/aggregator/mld_d_dec_map_bld.f90 similarity index 100% rename from mlprec/impl/mld_d_dec_map_bld.f90 rename to mlprec/impl/aggregator/mld_d_dec_map_bld.f90 diff --git a/mlprec/impl/aggregator/mld_d_hyb_map_bld.f90 b/mlprec/impl/aggregator/mld_d_hyb_map_bld.f90 new file mode 100644 index 00000000..a0f7f7cb --- /dev/null +++ b/mlprec/impl/aggregator/mld_d_hyb_map_bld.f90 @@ -0,0 +1,328 @@ +! +! +! 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. +! +! +! +! File: mld_d_hyb_map__bld.f90 +! +! Subroutine: mld_d_hyb_map_bld +! Version: real +! +! This routine builds the tentative prolongator based on 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. +! +! Note: upon exit +! +! Arguments: +! a - type(psb_dspmat_type). +! The sparse matrix structure containing the local part of the +! matrix to be preconditioned. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_dprec_type), input/output. +! The preconditioner data structure; upon exit it contains +! the multilevel hierarchy of prolongators, restrictors +! and coarse matrices. +! info - integer, output. +! Error code. +! +! +! +subroutine mld_d_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_hyb_map_bld + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:), tmpaggr(:) + real(psb_dpk_), allocatable :: val(:), diag(:) + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt + type(psb_d_csr_sparse_mat) :: acsr, muij, s_neigh + type(psb_d_coo_sparse_mat) :: s_neigh_coo + real(psb_dpk_) :: cpling, tcl + logical :: disjoint + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if (psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_hyb_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + nr = a%get_nrows() + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),icol(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + diag = a%get_diag(info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + ! + ! Phase zero: compute muij + ! + call a%cp_to(muij) + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) + end do + end do + !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! + ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 + ! + call s_neigh_coo%allocate(nr,nr,muij%get_nzeros()) + ip = 0 + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + ip = ip + 1 + s_neigh_coo%ia(ip) = i + s_neigh_coo%ja(ip) = j + if (real(muij%val(k)) >= theta) then + s_neigh_coo%val(ip) = done + else + s_neigh_coo%val(ip) = -done + end if + end do + end do + !write(*,*) 'S_NEIGH: ',nr,ip + call s_neigh_coo%set_nzeros(ip) + call s_neigh%mv_from_coo(s_neigh_coo,info) + + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = muij%irp(i+1) - muij%irp(i) + end do + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if + + + ! + ! Phase one: Start with disjoint groups. + ! + naggr = 0 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Get the 1-neighbourhood of I + ! + ip1 = s_neigh%irp(i) + nz = s_neigh%irp(i+1)-ip1 + ! + ! If the neighbourhood only contains I, skip it + ! + if (nz ==0) then + ilaggr(i) = 0 + cycle step1 + end if + if ((nz==1).and.(s_neigh%ja(ip1)==i)) then + ilaggr(i) = 0 + cycle step1 + end if + ! + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate. + ! + nzcnt = count(real(s_neigh%val(ip1:ip1+nz-1)) > 0) + icol(1:nzcnt) = pack(s_neigh%ja(ip1:ip1+nz-1),(real(s_neigh%val(ip1:ip1+nz-1)) > 0)) + disjoint = all(ilaggr(icol(1:nzcnt)) == -(nr+1)) + if (disjoint) then + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, nzcnt + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 + + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: join the neighbours + ! + tmpaggr = ilaggr + step2: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = dzero + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if ( (tmpaggr(j) > 0).and. (real(muij%val(k)) > cpling)& + & .and.(real(s_neigh%val(k))>0)) then + ip = k + cpling = muij%val(k) + end if + end if + enddo + if (ip > 0) then + ilaggr(i) = ilaggr(s_neigh%ja(ip)) + end if + end if + end do step2 + + + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) < 0) then + ! + ! Find its strongly connected neighbourhood not + ! already aggregated, and make it into a new aggregate. + ! + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if (ilaggr(j) < 0) then + ip = ip + 1 + icol(ip) = j + end if + end if + enddo + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do + end if + end if + end do step3 + + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_d_hyb_map_bld + diff --git a/mlprec/impl/aggregator/mld_d_hybrid_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_d_hybrid_aggregator_tprol.f90 new file mode 100644 index 00000000..43506b48 --- /dev/null +++ b/mlprec/impl/aggregator/mld_d_hybrid_aggregator_tprol.f90 @@ -0,0 +1,122 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_d_hybrid_aggregator_tprol.f90 +! +! Subroutine: mld_d_hybrid_aggregator_tprol +! Version: real +! +! +! This routine is mainly an interface to hyb_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_d_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_dspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_dspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_d_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_d_hybrid_aggregator_mod, mld_protect_name => mld_d_hybrid_aggregator_build_tprol + use mld_d_inner_mod + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_d_hybrid_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(p%parms%ml_cycle,'Multilevel cycle',& + & mld_mult_ml_,is_legal_ml_cycle) + call mld_check_def(p%parms%par_aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_par_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + call mld_hyb_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_d_hybrid_aggregator_build_tprol diff --git a/mlprec/impl/aggregator/mld_d_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_d_map_to_tprol.f90 new file mode 100644 index 00000000..1688b979 --- /dev/null +++ b/mlprec/impl/aggregator/mld_d_map_to_tprol.f90 @@ -0,0 +1,150 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_d_map_to_tprol.f90 +! +! Subroutine: mld_d_map_to_tprol +! Version: real +! +! This routine uses a mapping from the row indices of the fine-level matrix +! to the row indices of the coarse-level matrix to build a tentative +! prolongator, i.e. a piecewise constant operator. +! This is later used to build the final operator; the code has been refactored here +! to be shared among all the methods that provide the tentative prolongator +! through a simple integer mapping. +! +! The aggregation algorithm is a parallel version of that described in +! * M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! * P. Vanek, J. Mandel and M. Brezina, Algebraic Multigrid by Smoothed +! Aggregation for Second and Fourth Order Elliptic Problems, Computing, 56 +! (1996), 179-196. +! For more details see +! 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. +! +! +! Arguments: +! aggr_type - integer, input. +! The scalar used to identify the aggregation algorithm. +! theta - real, input. +! The aggregation threshold used in the aggregation algorithm. +! a - type(psb_dspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! ilaggr - integer, dimension(:), allocatable. +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable. +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_dspmat_type). +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_d_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_map_to_tprol + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr + type(psb_d_coo_sparse_mat) :: tmpcoo + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_map_to_tprol' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 + call psb_halo(ilaggr,desc_a,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') + goto 9999 + end if + + call tmpcoo%allocate(ncol,ntaggr,ncol) + do i=1,ncol + tmpcoo%val(i) = done + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) + end do + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_sorted() ! At this point this is in row-major + call op_prol%mv_from(tmpcoo) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_d_map_to_tprol diff --git a/mlprec/impl/aggregator/mld_d_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_d_symdec_aggregator_tprol.f90 new file mode 100644 index 00000000..8133b712 --- /dev/null +++ b/mlprec/impl/aggregator/mld_d_symdec_aggregator_tprol.f90 @@ -0,0 +1,141 @@ +! +! +! 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. +! +! +! File: mld_d_symdec_aggregator_tprol.f90 +! +! Subroutine: mld_d_symdec_aggregator_tprol +! Version: real +! +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! +! Arguments: +! p - type(mld_d_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_dspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_dspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_d_symdec_aggregator_mod, mld_protect_name => mld_d_symdec_aggregator_build_tprol + use mld_d_inner_mod + 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 + + ! Local variables + type(psb_dspmat_type) :: atmp, atrans + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr, nr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_d_symdec_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + + nr = a%get_nrows() + call a%csclip(atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atmp%transp(atrans) + if (info == psb_success_) call atrans%cscnv(info,type='COO') + if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atrans%free() + if (info == psb_success_) call atmp%cscnv(info,type='CSR') + + if (info == psb_success_) & + & call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call atmp%free() + + if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_d_symdec_aggregator_build_tprol diff --git a/mlprec/impl/mld_daggrmat_biz_asb.f90 b/mlprec/impl/aggregator/mld_daggrmat_biz_asb.f90 similarity index 100% rename from mlprec/impl/mld_daggrmat_biz_asb.f90 rename to mlprec/impl/aggregator/mld_daggrmat_biz_asb.f90 diff --git a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_daggrmat_minnrg_asb.f90 similarity index 100% rename from mlprec/impl/mld_daggrmat_minnrg_asb.f90 rename to mlprec/impl/aggregator/mld_daggrmat_minnrg_asb.f90 diff --git a/mlprec/impl/mld_daggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_daggrmat_nosmth_asb.f90 similarity index 100% rename from mlprec/impl/mld_daggrmat_nosmth_asb.f90 rename to mlprec/impl/aggregator/mld_daggrmat_nosmth_asb.f90 diff --git a/mlprec/impl/mld_daggrmat_smth_asb.f90 b/mlprec/impl/aggregator/mld_daggrmat_smth_asb.f90 similarity index 100% rename from mlprec/impl/mld_daggrmat_smth_asb.f90 rename to mlprec/impl/aggregator/mld_daggrmat_smth_asb.f90 diff --git a/mlprec/impl/aggregator/mld_s_base_aggregator_mat_asb.f90 b/mlprec/impl/aggregator/mld_s_base_aggregator_mat_asb.f90 new file mode 100644 index 00000000..4464cc84 --- /dev/null +++ b/mlprec/impl/aggregator/mld_s_base_aggregator_mat_asb.f90 @@ -0,0 +1,225 @@ +! +! +! 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. +! +! +! File: mld_s_base_aggregator_mat_asb.f90 +! +! Subroutine: mld_s_base_aggregator_mat_asb +! Version: real +! +! This routine builds the matrix associated to the current level of the +! multilevel preconditioner from the matrix associated to the previous level, +! by using the user-specified aggregation technique (therefore, it also builds the +! prolongation and restriction operators mapping the current level to the +! previous one and vice versa). +! The current level is regarded as the coarse one, while the previous as +! the fine one. This is in agreement with the fact that the routine is called, +! by mld_mlprec_bld, only on levels >=2. +! The coarse-level matrix A_C is built from a fine-level matrix A +! by using the Galerkin approach, i.e. +! +! A_C = P_C^T A P_C, +! +! where P_C is a prolongator from the coarse level to the fine one. +! +! A mapping from the nodes of the adjacency graph of A to the nodes of the +! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine. +! The prolongator P_C is built here from this mapping, according to the +! value of p%iprcparm(mld_aggr_kind_), specified by the user through +! mld_sprecinit and mld_zprecset. +! On output from this routine the entries of AC, op_prol, op_restr +! are still in "global numbering" mode; this is fixed in the calling routine +! mld_s_lev_aggrmat_asb. +! +! Currently four different prolongators are implemented, corresponding to +! four aggregation algorithms: +! 1. un-smoothed aggregation, +! 2. smoothed aggregation, +! 3. "bizarre" aggregation. +! 4. minimum energy +! 1. The non-smoothed aggregation uses as prolongator the piecewise constant +! interpolation operator corresponding to the fine-to-coarse level mapping built +! by p%aggr%bld_tprol. This is called tentative prolongator. +! 2. The smoothed aggregation uses as prolongator the operator obtained by applying +! a damped Jacobi smoother to the tentative prolongator. +! 3. The "bizarre" aggregation uses a prolongator proposed by the authors of MLD2P4. +! This prolongator still requires a deep analysis and testing and its use is +! not recommended. +! 4. Minimum energy aggregation +! +! For more details see +! 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. +! M. Sala, R. Tuminaro: A new Petrov-Galerkin smoothed aggregation preconditioner +! for nonsymmetric linear systems, SIAM J. Sci. Comput., 31(1):143-166 (2008) +! +! +! The main structure is: +! 1. Perform sanity checks; +! 2. Compute prolongator/restrictor/AC +! +! +! Arguments: +! ag - type(mld_s_base_aggregator_type), input/output. +! The aggregator object +! parms - type(mld_sml_parms), input +! The aggregation parameters +! a - type(psb_sspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! The 'one-level' data structure that will contain the local +! part of the matrix to be built as well as the information +! concerning the prolongator and its transpose. +! ilaggr - integer, dimension(:), input +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that the indices +! are assumed to be shifted so as to make sure the ranges on +! the various processes do not overlap. +! nlaggr - integer, dimension(:) input +! nlaggr(i) contains the aggregates held by process i. +! ac - type(psb_sspmat_type), output +! The coarse matrix on output +! +! op_prol - type(psb_sspmat_type), input/output +! The tentative prolongator on input, the computed prolongator on output +! +! op_restr - type(psb_sspmat_type), output +! The restrictor operator; normally, it is the transpose of the prolongator. +! +! info - integer, output. +! Error code. +! +subroutine mld_s_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info) + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_base_aggregator_mat_asb + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + type(psb_s_coo_sparse_mat) :: acoo, bcoo + type(psb_s_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl,ntaggr + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_s_base_aggregator_mat_asb' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%aggr_kind,'Smoother',& + & mld_smooth_prol_,is_legal_ml_aggr_kind) + call mld_check_def(parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + call mld_check_def(parms%aggr_filter,'Use filtered matrix',& + & mld_no_filter_mat_,is_legal_aggr_filter) + call mld_check_def(parms%smoother_pos,'smooth_pos',& + & mld_pre_smooth_,is_legal_ml_smooth_pos) + call mld_check_def(parms%aggr_omega_alg,'Omega Alg.',& + & mld_eig_est_,is_legal_ml_aggr_omega_alg) + call mld_check_def(parms%aggr_eig,'Eigenvalue estimate',& + & mld_max_norm_,is_legal_ml_aggr_eig) + call mld_check_def(parms%aggr_omega_val,'Omega',szero,is_legal_s_omega) + + ! + ! Build the coarse-level matrix from the fine-level one, starting from + ! the mapping defined by mld_aggrmap_bld and applying the aggregation + ! algorithm specified by p%iprcparm(mld_aggr_kind_) + ! + select case (parms%aggr_kind) + case (mld_no_smooth_) + + call mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & parms,ac,op_prol,op_restr,info) + + case(mld_smooth_prol_) + + call mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_biz_prol_) + + call mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + +end subroutine mld_s_base_aggregator_mat_asb diff --git a/mlprec/impl/aggregator/mld_s_base_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_s_base_aggregator_tprol.f90 new file mode 100644 index 00000000..bf83a553 --- /dev/null +++ b/mlprec/impl/aggregator/mld_s_base_aggregator_tprol.f90 @@ -0,0 +1,124 @@ +! +! +! 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. +! +! +! File: mld_s_base_aggregator_tprol.f90 +! +! Subroutine: mld_s_base_aggregator_tprol +! Version: real +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_s_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_sspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_sspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_s_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod +! use mld_s_base_aggregator_mod + use mld_s_inner_mod, mld_protect_name => mld_s_base_aggregator_build_tprol + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_s_base_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) + + + call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_s_base_aggregator_build_tprol diff --git a/mlprec/impl/mld_s_dec_map_bld.f90 b/mlprec/impl/aggregator/mld_s_dec_map_bld.f90 similarity index 100% rename from mlprec/impl/mld_s_dec_map_bld.f90 rename to mlprec/impl/aggregator/mld_s_dec_map_bld.f90 diff --git a/mlprec/impl/aggregator/mld_s_hyb_map_bld.f90 b/mlprec/impl/aggregator/mld_s_hyb_map_bld.f90 new file mode 100644 index 00000000..8331ae05 --- /dev/null +++ b/mlprec/impl/aggregator/mld_s_hyb_map_bld.f90 @@ -0,0 +1,328 @@ +! +! +! 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. +! +! +! +! File: mld_s_hyb_map__bld.f90 +! +! Subroutine: mld_s_hyb_map_bld +! Version: real +! +! This routine builds the tentative prolongator based on 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. +! +! Note: upon exit +! +! Arguments: +! a - type(psb_sspmat_type). +! The sparse matrix structure containing the local part of the +! matrix to be preconditioned. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_sprec_type), input/output. +! The preconditioner data structure; upon exit it contains +! the multilevel hierarchy of prolongators, restrictors +! and coarse matrices. +! info - integer, output. +! Error code. +! +! +! +subroutine mld_s_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_hyb_map_bld + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:), tmpaggr(:) + real(psb_spk_), allocatable :: val(:), diag(:) + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt + type(psb_s_csr_sparse_mat) :: acsr, muij, s_neigh + type(psb_s_coo_sparse_mat) :: s_neigh_coo + real(psb_spk_) :: cpling, tcl + logical :: disjoint + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if (psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_hyb_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + nr = a%get_nrows() + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),icol(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + diag = a%get_diag(info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + ! + ! Phase zero: compute muij + ! + call a%cp_to(muij) + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) + end do + end do + !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! + ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 + ! + call s_neigh_coo%allocate(nr,nr,muij%get_nzeros()) + ip = 0 + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + ip = ip + 1 + s_neigh_coo%ia(ip) = i + s_neigh_coo%ja(ip) = j + if (real(muij%val(k)) >= theta) then + s_neigh_coo%val(ip) = sone + else + s_neigh_coo%val(ip) = -sone + end if + end do + end do + !write(*,*) 'S_NEIGH: ',nr,ip + call s_neigh_coo%set_nzeros(ip) + call s_neigh%mv_from_coo(s_neigh_coo,info) + + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = muij%irp(i+1) - muij%irp(i) + end do + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if + + + ! + ! Phase one: Start with disjoint groups. + ! + naggr = 0 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Get the 1-neighbourhood of I + ! + ip1 = s_neigh%irp(i) + nz = s_neigh%irp(i+1)-ip1 + ! + ! If the neighbourhood only contains I, skip it + ! + if (nz ==0) then + ilaggr(i) = 0 + cycle step1 + end if + if ((nz==1).and.(s_neigh%ja(ip1)==i)) then + ilaggr(i) = 0 + cycle step1 + end if + ! + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate. + ! + nzcnt = count(real(s_neigh%val(ip1:ip1+nz-1)) > 0) + icol(1:nzcnt) = pack(s_neigh%ja(ip1:ip1+nz-1),(real(s_neigh%val(ip1:ip1+nz-1)) > 0)) + disjoint = all(ilaggr(icol(1:nzcnt)) == -(nr+1)) + if (disjoint) then + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, nzcnt + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 + + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: join the neighbours + ! + tmpaggr = ilaggr + step2: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = szero + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if ( (tmpaggr(j) > 0).and. (real(muij%val(k)) > cpling)& + & .and.(real(s_neigh%val(k))>0)) then + ip = k + cpling = muij%val(k) + end if + end if + enddo + if (ip > 0) then + ilaggr(i) = ilaggr(s_neigh%ja(ip)) + end if + end if + end do step2 + + + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) < 0) then + ! + ! Find its strongly connected neighbourhood not + ! already aggregated, and make it into a new aggregate. + ! + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if (ilaggr(j) < 0) then + ip = ip + 1 + icol(ip) = j + end if + end if + enddo + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do + end if + end if + end do step3 + + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_hyb_map_bld + diff --git a/mlprec/impl/aggregator/mld_s_hybrid_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_s_hybrid_aggregator_tprol.f90 new file mode 100644 index 00000000..806f3184 --- /dev/null +++ b/mlprec/impl/aggregator/mld_s_hybrid_aggregator_tprol.f90 @@ -0,0 +1,122 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_s_hybrid_aggregator_tprol.f90 +! +! Subroutine: mld_s_hybrid_aggregator_tprol +! Version: real +! +! +! This routine is mainly an interface to hyb_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_s_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_sspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_sspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_s_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_s_hybrid_aggregator_mod, mld_protect_name => mld_s_hybrid_aggregator_build_tprol + use mld_s_inner_mod + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_s_hybrid_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(p%parms%ml_cycle,'Multilevel cycle',& + & mld_mult_ml_,is_legal_ml_cycle) + call mld_check_def(p%parms%par_aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_par_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) + call mld_hyb_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_s_hybrid_aggregator_build_tprol diff --git a/mlprec/impl/aggregator/mld_s_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_s_map_to_tprol.f90 new file mode 100644 index 00000000..0b41787c --- /dev/null +++ b/mlprec/impl/aggregator/mld_s_map_to_tprol.f90 @@ -0,0 +1,150 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_s_map_to_tprol.f90 +! +! Subroutine: mld_s_map_to_tprol +! Version: real +! +! This routine uses a mapping from the row indices of the fine-level matrix +! to the row indices of the coarse-level matrix to build a tentative +! prolongator, i.e. a piecewise constant operator. +! This is later used to build the final operator; the code has been refactored here +! to be shared among all the methods that provide the tentative prolongator +! through a simple integer mapping. +! +! The aggregation algorithm is a parallel version of that described in +! * M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! * P. Vanek, J. Mandel and M. Brezina, Algebraic Multigrid by Smoothed +! Aggregation for Second and Fourth Order Elliptic Problems, Computing, 56 +! (1996), 179-196. +! For more details see +! 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. +! +! +! Arguments: +! aggr_type - integer, input. +! The scalar used to identify the aggregation algorithm. +! theta - real, input. +! The aggregation threshold used in the aggregation algorithm. +! a - type(psb_sspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! ilaggr - integer, dimension(:), allocatable. +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable. +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_sspmat_type). +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_s_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_map_to_tprol + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr + type(psb_s_coo_sparse_mat) :: tmpcoo + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_map_to_tprol' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 + call psb_halo(ilaggr,desc_a,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') + goto 9999 + end if + + call tmpcoo%allocate(ncol,ntaggr,ncol) + do i=1,ncol + tmpcoo%val(i) = sone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) + end do + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_sorted() ! At this point this is in row-major + call op_prol%mv_from(tmpcoo) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_map_to_tprol diff --git a/mlprec/impl/aggregator/mld_s_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_s_symdec_aggregator_tprol.f90 new file mode 100644 index 00000000..b73fe100 --- /dev/null +++ b/mlprec/impl/aggregator/mld_s_symdec_aggregator_tprol.f90 @@ -0,0 +1,141 @@ +! +! +! 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. +! +! +! File: mld_s_symdec_aggregator_tprol.f90 +! +! Subroutine: mld_s_symdec_aggregator_tprol +! Version: real +! +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! +! Arguments: +! p - type(mld_s_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_sspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_sspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_s_symdec_aggregator_mod, mld_protect_name => mld_s_symdec_aggregator_build_tprol + use mld_s_inner_mod + 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 + + ! Local variables + type(psb_sspmat_type) :: atmp, atrans + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr, nr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_s_symdec_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) + + nr = a%get_nrows() + call a%csclip(atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atmp%transp(atrans) + if (info == psb_success_) call atrans%cscnv(info,type='COO') + if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atrans%free() + if (info == psb_success_) call atmp%cscnv(info,type='CSR') + + if (info == psb_success_) & + & call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call atmp%free() + + if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_s_symdec_aggregator_build_tprol diff --git a/mlprec/impl/mld_saggrmat_biz_asb.f90 b/mlprec/impl/aggregator/mld_saggrmat_biz_asb.f90 similarity index 100% rename from mlprec/impl/mld_saggrmat_biz_asb.f90 rename to mlprec/impl/aggregator/mld_saggrmat_biz_asb.f90 diff --git a/mlprec/impl/mld_saggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_saggrmat_minnrg_asb.f90 similarity index 100% rename from mlprec/impl/mld_saggrmat_minnrg_asb.f90 rename to mlprec/impl/aggregator/mld_saggrmat_minnrg_asb.f90 diff --git a/mlprec/impl/mld_saggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_saggrmat_nosmth_asb.f90 similarity index 100% rename from mlprec/impl/mld_saggrmat_nosmth_asb.f90 rename to mlprec/impl/aggregator/mld_saggrmat_nosmth_asb.f90 diff --git a/mlprec/impl/mld_saggrmat_smth_asb.f90 b/mlprec/impl/aggregator/mld_saggrmat_smth_asb.f90 similarity index 100% rename from mlprec/impl/mld_saggrmat_smth_asb.f90 rename to mlprec/impl/aggregator/mld_saggrmat_smth_asb.f90 diff --git a/mlprec/impl/aggregator/mld_z_base_aggregator_mat_asb.f90 b/mlprec/impl/aggregator/mld_z_base_aggregator_mat_asb.f90 new file mode 100644 index 00000000..5bda2a48 --- /dev/null +++ b/mlprec/impl/aggregator/mld_z_base_aggregator_mat_asb.f90 @@ -0,0 +1,225 @@ +! +! +! 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. +! +! +! File: mld_z_base_aggregator_mat_asb.f90 +! +! Subroutine: mld_z_base_aggregator_mat_asb +! Version: complex +! +! This routine builds the matrix associated to the current level of the +! multilevel preconditioner from the matrix associated to the previous level, +! by using the user-specified aggregation technique (therefore, it also builds the +! prolongation and restriction operators mapping the current level to the +! previous one and vice versa). +! The current level is regarded as the coarse one, while the previous as +! the fine one. This is in agreement with the fact that the routine is called, +! by mld_mlprec_bld, only on levels >=2. +! The coarse-level matrix A_C is built from a fine-level matrix A +! by using the Galerkin approach, i.e. +! +! A_C = P_C^T A P_C, +! +! where P_C is a prolongator from the coarse level to the fine one. +! +! A mapping from the nodes of the adjacency graph of A to the nodes of the +! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine. +! The prolongator P_C is built here from this mapping, according to the +! value of p%iprcparm(mld_aggr_kind_), specified by the user through +! mld_zprecinit and mld_zprecset. +! On output from this routine the entries of AC, op_prol, op_restr +! are still in "global numbering" mode; this is fixed in the calling routine +! mld_z_lev_aggrmat_asb. +! +! Currently four different prolongators are implemented, corresponding to +! four aggregation algorithms: +! 1. un-smoothed aggregation, +! 2. smoothed aggregation, +! 3. "bizarre" aggregation. +! 4. minimum energy +! 1. The non-smoothed aggregation uses as prolongator the piecewise constant +! interpolation operator corresponding to the fine-to-coarse level mapping built +! by p%aggr%bld_tprol. This is called tentative prolongator. +! 2. The smoothed aggregation uses as prolongator the operator obtained by applying +! a damped Jacobi smoother to the tentative prolongator. +! 3. The "bizarre" aggregation uses a prolongator proposed by the authors of MLD2P4. +! This prolongator still requires a deep analysis and testing and its use is +! not recommended. +! 4. Minimum energy aggregation +! +! For more details see +! 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. +! M. Sala, R. Tuminaro: A new Petrov-Galerkin smoothed aggregation preconditioner +! for nonsymmetric linear systems, SIAM J. Sci. Comput., 31(1):143-166 (2008) +! +! +! The main structure is: +! 1. Perform sanity checks; +! 2. Compute prolongator/restrictor/AC +! +! +! Arguments: +! ag - type(mld_z_base_aggregator_type), input/output. +! The aggregator object +! parms - type(mld_dml_parms), input +! The aggregation parameters +! a - type(psb_zspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! The 'one-level' data structure that will contain the local +! part of the matrix to be built as well as the information +! concerning the prolongator and its transpose. +! ilaggr - integer, dimension(:), input +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that the indices +! are assumed to be shifted so as to make sure the ranges on +! the various processes do not overlap. +! nlaggr - integer, dimension(:) input +! nlaggr(i) contains the aggregates held by process i. +! ac - type(psb_zspmat_type), output +! The coarse matrix on output +! +! op_prol - type(psb_zspmat_type), input/output +! The tentative prolongator on input, the computed prolongator on output +! +! op_restr - type(psb_zspmat_type), output +! The restrictor operator; normally, it is the transpose of the prolongator. +! +! info - integer, output. +! Error code. +! +subroutine mld_z_base_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,op_prol,op_restr,info) + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_base_aggregator_mat_asb + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + type(psb_z_coo_sparse_mat) :: acoo, bcoo + type(psb_z_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl,ntaggr + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_z_base_aggregator_mat_asb' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%aggr_kind,'Smoother',& + & mld_smooth_prol_,is_legal_ml_aggr_kind) + call mld_check_def(parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + call mld_check_def(parms%aggr_filter,'Use filtered matrix',& + & mld_no_filter_mat_,is_legal_aggr_filter) + call mld_check_def(parms%smoother_pos,'smooth_pos',& + & mld_pre_smooth_,is_legal_ml_smooth_pos) + call mld_check_def(parms%aggr_omega_alg,'Omega Alg.',& + & mld_eig_est_,is_legal_ml_aggr_omega_alg) + call mld_check_def(parms%aggr_eig,'Eigenvalue estimate',& + & mld_max_norm_,is_legal_ml_aggr_eig) + call mld_check_def(parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega) + + ! + ! Build the coarse-level matrix from the fine-level one, starting from + ! the mapping defined by mld_aggrmap_bld and applying the aggregation + ! algorithm specified by p%iprcparm(mld_aggr_kind_) + ! + select case (parms%aggr_kind) + case (mld_no_smooth_) + + call mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & parms,ac,op_prol,op_restr,info) + + case(mld_smooth_prol_) + + call mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_biz_prol_) + + call mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & + & parms,ac,op_prol,op_restr,info) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + +end subroutine mld_z_base_aggregator_mat_asb diff --git a/mlprec/impl/aggregator/mld_z_base_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_z_base_aggregator_tprol.f90 new file mode 100644 index 00000000..81d8910f --- /dev/null +++ b/mlprec/impl/aggregator/mld_z_base_aggregator_tprol.f90 @@ -0,0 +1,124 @@ +! +! +! 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. +! +! +! File: mld_z_base_aggregator_tprol.f90 +! +! Subroutine: mld_z_base_aggregator_tprol +! Version: complex +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_z_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_zspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_zspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_z_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod +! use mld_z_base_aggregator_mod + use mld_z_inner_mod, mld_protect_name => mld_z_base_aggregator_build_tprol + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_z_base_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + + + call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_z_base_aggregator_build_tprol diff --git a/mlprec/impl/mld_z_dec_map_bld.f90 b/mlprec/impl/aggregator/mld_z_dec_map_bld.f90 similarity index 100% rename from mlprec/impl/mld_z_dec_map_bld.f90 rename to mlprec/impl/aggregator/mld_z_dec_map_bld.f90 diff --git a/mlprec/impl/aggregator/mld_z_hyb_map_bld.f90 b/mlprec/impl/aggregator/mld_z_hyb_map_bld.f90 new file mode 100644 index 00000000..cc795465 --- /dev/null +++ b/mlprec/impl/aggregator/mld_z_hyb_map_bld.f90 @@ -0,0 +1,328 @@ +! +! +! 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. +! +! +! +! File: mld_z_hyb_map__bld.f90 +! +! Subroutine: mld_z_hyb_map_bld +! Version: complex +! +! This routine builds the tentative prolongator based on 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. +! +! Note: upon exit +! +! Arguments: +! a - type(psb_zspmat_type). +! The sparse matrix structure containing the local part of the +! matrix to be preconditioned. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_zprec_type), input/output. +! The preconditioner data structure; upon exit it contains +! the multilevel hierarchy of prolongators, restrictors +! and coarse matrices. +! info - integer, output. +! Error code. +! +! +! +subroutine mld_z_hyb_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_hyb_map_bld + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:), tmpaggr(:) + complex(psb_dpk_), allocatable :: val(:), diag(:) + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip, ip1,nzcnt + type(psb_z_csr_sparse_mat) :: acsr, muij, s_neigh + type(psb_z_coo_sparse_mat) :: s_neigh_coo + real(psb_dpk_) :: cpling, tcl + logical :: disjoint + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if (psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_hyb_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + nr = a%get_nrows() + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),icol(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + diag = a%get_diag(info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + ! + ! Phase zero: compute muij + ! + call a%cp_to(muij) + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) + end do + end do + !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! + ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 + ! + call s_neigh_coo%allocate(nr,nr,muij%get_nzeros()) + ip = 0 + do i=1, nr + do k=muij%irp(i),muij%irp(i+1)-1 + j = muij%ja(k) + ip = ip + 1 + s_neigh_coo%ia(ip) = i + s_neigh_coo%ja(ip) = j + if (real(muij%val(k)) >= theta) then + s_neigh_coo%val(ip) = done + else + s_neigh_coo%val(ip) = -done + end if + end do + end do + !write(*,*) 'S_NEIGH: ',nr,ip + call s_neigh_coo%set_nzeros(ip) + call s_neigh%mv_from_coo(s_neigh_coo,info) + + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = muij%irp(i+1) - muij%irp(i) + end do + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if + + + ! + ! Phase one: Start with disjoint groups. + ! + naggr = 0 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Get the 1-neighbourhood of I + ! + ip1 = s_neigh%irp(i) + nz = s_neigh%irp(i+1)-ip1 + ! + ! If the neighbourhood only contains I, skip it + ! + if (nz ==0) then + ilaggr(i) = 0 + cycle step1 + end if + if ((nz==1).and.(s_neigh%ja(ip1)==i)) then + ilaggr(i) = 0 + cycle step1 + end if + ! + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate. + ! + nzcnt = count(real(s_neigh%val(ip1:ip1+nz-1)) > 0) + icol(1:nzcnt) = pack(s_neigh%ja(ip1:ip1+nz-1),(real(s_neigh%val(ip1:ip1+nz-1)) > 0)) + disjoint = all(ilaggr(icol(1:nzcnt)) == -(nr+1)) + if (disjoint) then + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, nzcnt + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 + + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: join the neighbours + ! + tmpaggr = ilaggr + step2: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = dzero + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if ( (tmpaggr(j) > 0).and. (real(muij%val(k)) > cpling)& + & .and.(real(s_neigh%val(k))>0)) then + ip = k + cpling = muij%val(k) + end if + end if + enddo + if (ip > 0) then + ilaggr(i) = ilaggr(s_neigh%ja(ip)) + end if + end if + end do step2 + + + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) < 0) then + ! + ! Find its strongly connected neighbourhood not + ! already aggregated, and make it into a new aggregate. + ! + ip = 0 + do k=s_neigh%irp(i), s_neigh%irp(i+1)-1 + j = s_neigh%ja(k) + if ((1<=j).and.(j<=nr)) then + if (ilaggr(j) < 0) then + ip = ip + 1 + icol(ip) = j + end if + end if + enddo + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do + end if + end if + end do step3 + + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_hyb_map_bld + diff --git a/mlprec/impl/aggregator/mld_z_hybrid_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_z_hybrid_aggregator_tprol.f90 new file mode 100644 index 00000000..d0a83ea4 --- /dev/null +++ b/mlprec/impl/aggregator/mld_z_hybrid_aggregator_tprol.f90 @@ -0,0 +1,122 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_z_hybrid_aggregator_tprol.f90 +! +! Subroutine: mld_z_hybrid_aggregator_tprol +! Version: complex +! +! +! This routine is mainly an interface to hyb_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! Arguments: +! p - type(mld_z_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_zspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_zspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_z_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_z_hybrid_aggregator_mod, mld_protect_name => mld_z_hybrid_aggregator_build_tprol + use mld_z_inner_mod + 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 + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_z_hybrid_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(p%parms%ml_cycle,'Multilevel cycle',& + & mld_mult_ml_,is_legal_ml_cycle) + call mld_check_def(p%parms%par_aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_par_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + call mld_hyb_map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + + call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_z_hybrid_aggregator_build_tprol diff --git a/mlprec/impl/aggregator/mld_z_map_to_tprol.f90 b/mlprec/impl/aggregator/mld_z_map_to_tprol.f90 new file mode 100644 index 00000000..5a9491e9 --- /dev/null +++ b/mlprec/impl/aggregator/mld_z_map_to_tprol.f90 @@ -0,0 +1,150 @@ +! +! +! 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 +! +! Salvatore Filippone Cranfield University, UK +! 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. +! +! +! File: mld_z_map_to_tprol.f90 +! +! Subroutine: mld_z_map_to_tprol +! Version: complex +! +! This routine uses a mapping from the row indices of the fine-level matrix +! to the row indices of the coarse-level matrix to build a tentative +! prolongator, i.e. a piecewise constant operator. +! This is later used to build the final operator; the code has been refactored here +! to be shared among all the methods that provide the tentative prolongator +! through a simple integer mapping. +! +! The aggregation algorithm is a parallel version of that described in +! * M. Brezina and P. Vanek, A black-box iterative solver based on a +! two-level Schwarz method, Computing, 63 (1999), 233-263. +! * P. Vanek, J. Mandel and M. Brezina, Algebraic Multigrid by Smoothed +! Aggregation for Second and Fourth Order Elliptic Problems, Computing, 56 +! (1996), 179-196. +! For more details see +! 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. +! +! +! Arguments: +! aggr_type - integer, input. +! The scalar used to identify the aggregation algorithm. +! theta - real, input. +! The aggregation threshold used in the aggregation algorithm. +! a - type(psb_zspmat_type), input. +! The sparse matrix structure containing the local part of +! the fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the fine-level matrix. +! ilaggr - integer, dimension(:), allocatable. +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable. +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_zspmat_type). +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_z_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_map_to_tprol + + implicit none + + ! Arguments + 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 + + ! Local variables + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr + type(psb_z_coo_sparse_mat) :: tmpcoo + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_map_to_tprol' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 + call psb_halo(ilaggr,desc_a,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') + goto 9999 + end if + + call tmpcoo%allocate(ncol,ntaggr,ncol) + do i=1,ncol + tmpcoo%val(i) = zone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) + end do + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_sorted() ! At this point this is in row-major + call op_prol%mv_from(tmpcoo) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_map_to_tprol diff --git a/mlprec/impl/aggregator/mld_z_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_z_symdec_aggregator_tprol.f90 new file mode 100644 index 00000000..7a8d5967 --- /dev/null +++ b/mlprec/impl/aggregator/mld_z_symdec_aggregator_tprol.f90 @@ -0,0 +1,141 @@ +! +! +! 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. +! +! +! File: mld_z_symdec_aggregator_tprol.f90 +! +! Subroutine: mld_z_symdec_aggregator_tprol +! Version: complex +! +! +! This routine is mainly an interface to dec_map_bld where the real work is performed. +! It takes care of some consistency checking, and calls map_to_tprol, which is +! refactored and shared among all the aggregation methods that produce a simple +! integer mapping. +! +! +! +! Arguments: +! p - type(mld_z_onelev_type), input/output. +! The 'one-level' data structure containing the control +! parameters and (eventually) coarse matrix and prolongator/restrictors. +! +! a - type(psb_zspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! ilaggr - integer, dimension(:), allocatable, output +! The mapping between the row indices of the coarse-level +! matrix and the row indices of the fine-level matrix. +! ilaggr(i)=j means that node i in the adjacency graph +! of the fine-level matrix is mapped onto node j in the +! adjacency graph of the coarse-level matrix. Note that on exit the indices +! will be shifted so as to make sure the ranges on the various processes do not +! overlap. +! nlaggr - integer, dimension(:), allocatable, output +! nlaggr(i) contains the aggregates held by process i. +! op_prol - type(psb_zspmat_type), output +! The tentative prolongator, based on ilaggr. +! +! info - integer, output. +! Error code. +! +subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_z_symdec_aggregator_mod, mld_protect_name => mld_z_symdec_aggregator_build_tprol + use mld_z_inner_mod + 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 + + ! Local variables + type(psb_zspmat_type) :: atmp, atrans + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ntaggr, nr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_z_symdec_aggregator_tprol' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + + nr = a%get_nrows() + call a%csclip(atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atmp%transp(atrans) + if (info == psb_success_) call atrans%cscnv(info,type='COO') + if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + call atmp%set_nrows(nr) + call atmp%set_ncols(nr) + if (info == psb_success_) call atrans%free() + if (info == psb_success_) call atmp%cscnv(info,type='CSR') + + if (info == psb_success_) & + & call mld_dec_map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call atmp%free() + + if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_z_symdec_aggregator_build_tprol diff --git a/mlprec/impl/mld_zaggrmat_biz_asb.f90 b/mlprec/impl/aggregator/mld_zaggrmat_biz_asb.f90 similarity index 100% rename from mlprec/impl/mld_zaggrmat_biz_asb.f90 rename to mlprec/impl/aggregator/mld_zaggrmat_biz_asb.f90 diff --git a/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 b/mlprec/impl/aggregator/mld_zaggrmat_minnrg_asb.f90 similarity index 100% rename from mlprec/impl/mld_zaggrmat_minnrg_asb.f90 rename to mlprec/impl/aggregator/mld_zaggrmat_minnrg_asb.f90 diff --git a/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 b/mlprec/impl/aggregator/mld_zaggrmat_nosmth_asb.f90 similarity index 100% rename from mlprec/impl/mld_zaggrmat_nosmth_asb.f90 rename to mlprec/impl/aggregator/mld_zaggrmat_nosmth_asb.f90 diff --git a/mlprec/impl/mld_zaggrmat_smth_asb.f90 b/mlprec/impl/aggregator/mld_zaggrmat_smth_asb.f90 similarity index 100% rename from mlprec/impl/mld_zaggrmat_smth_asb.f90 rename to mlprec/impl/aggregator/mld_zaggrmat_smth_asb.f90