From 9dab2a8c7c1f061fac2a9384b48b66e00796f9d8 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 20 Sep 2021 10:59:44 -0400 Subject: [PATCH] Further build new matching --- amgprec/amg_d_newmatch_aggregator_mod.F90 | 34 +++ amgprec/impl/aggregator/Makefile | 4 +- ...mg_d_newmatch_aggregator_inner_mat_asb.F90 | 166 +++++++++++++ .../amg_d_newmatch_aggregator_tprol.F90 | 20 +- .../amg_d_newmatch_spmm_bld_inner.F90 | 218 ++++++++++++++++++ .../aggregator/amg_d_newmatch_spmm_bld_ov.F90 | 169 ++++++++++++++ 6 files changed, 600 insertions(+), 11 deletions(-) create mode 100644 amgprec/impl/aggregator/amg_d_newmatch_aggregator_inner_mat_asb.F90 create mode 100644 amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_inner.F90 create mode 100644 amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_ov.F90 diff --git a/amgprec/amg_d_newmatch_aggregator_mod.F90 b/amgprec/amg_d_newmatch_aggregator_mod.F90 index 7266e09b..c55c5eb3 100644 --- a/amgprec/amg_d_newmatch_aggregator_mod.F90 +++ b/amgprec/amg_d_newmatch_aggregator_mod.F90 @@ -228,6 +228,40 @@ module amg_d_newmatch_aggregator_mod !!$ end subroutine amg_d_newmatch_unsmth_spmm_bld !!$ end interface + interface + subroutine amg_d_newmatch_spmm_bld_ov(a,desc_a,ilaggr,nlaggr,parms,& + & ac,desc_ac,op_prol,op_restr,t_prol,info) + import :: amg_d_newmatch_aggregator_type, psb_desc_type, psb_dspmat_type,& + & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data + implicit none + type(psb_dspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_ldspmat_type), intent(inout) :: t_prol + type(psb_dspmat_type), intent(inout) :: op_prol,ac, op_restr + type(psb_desc_type), intent(out) :: desc_ac + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_newmatch_spmm_bld_ov + end interface + + interface + subroutine amg_d_newmatch_spmm_bld_inner(a,desc_a,ilaggr,nlaggr,parms,& + & ac,desc_ac,op_prol,op_restr,t_prol,info) + import :: amg_d_newmatch_aggregator_type, psb_desc_type, psb_dspmat_type,& + & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data,& + & psb_d_csr_sparse_mat, psb_ld_csr_sparse_mat + implicit none + type(psb_d_csr_sparse_mat), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_ldspmat_type), intent(inout) :: t_prol + type(psb_dspmat_type), intent(inout) :: op_prol,ac, op_restr + type(psb_desc_type), intent(out) :: desc_ac + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_newmatch_spmm_bld_inner + end interface private :: is_legal_malg, is_legal_csize, is_legal_nsweeps, is_legal_nlevels diff --git a/amgprec/impl/aggregator/Makefile b/amgprec/impl/aggregator/Makefile index 7b664c48..49f76ed3 100644 --- a/amgprec/impl/aggregator/Makefile +++ b/amgprec/impl/aggregator/Makefile @@ -64,7 +64,9 @@ amg_d_newmatch_aggregator_mat_asb.o \ amg_d_newmatch_aggregator_inner_mat_asb.o \ amg_d_newmatch_aggregator_mat_bld.o \ amg_d_newmatch_aggregator_tprol.o \ -amg_d_newmatch_map_to_tprol.o +amg_d_newmatch_map_to_tprol.o \ +amg_d_newmatch_spmm_bld_inner.o \ +amg_d_newmatch_spmm_bld_ov.o MPCOBJS= #MatchBoxPC.o \ diff --git a/amgprec/impl/aggregator/amg_d_newmatch_aggregator_inner_mat_asb.F90 b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_inner_mat_asb.F90 new file mode 100644 index 00000000..e9a8b5de --- /dev/null +++ b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_inner_mat_asb.F90 @@ -0,0 +1,166 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! File: amg_d_newmatch_aggregator_mat_asb.f90 +! +! Subroutine: amg_d_newmatch_aggregator_mat_asb +! Version: real +! +! +! From a given AC to final format, generating DESC_AC. +! This is quite involved, because in the context of aggregation based +! on parallel matching we are building the matrix hierarchy within BLD_TPROL +! as we go, especially if we have multiple sweeps, hence this code is called +! in two completely different contexts: +! 1. Within bld_tprol for the internal hierarchy +! 2. Outside, from amg_hierarchy_bld +! The solution we have found is for bld_tprol to copy its output +! into special components ag%ac ag%desc_ac etc so that: +! 1. if they are allocated, it means that bld_tprol has been already invoked, we are in +! amg_hierarchy_bld and we only need to copy them +! 2. If they are not allocated, we are within bld_tprol, and we need to actually +! perform the various needed steps. +! +! Arguments: +! ag - type(amg_d_newmatch_aggregator_type), input/output. +! The aggregator object +! parms - type(amg_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), inout +! The coarse matrix +! desc_ac - type(psb_desc_type), output. +! 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. +! +! op_prol - type(psb_dspmat_type), input/output +! The tentative prolongator on input, the computed prolongator on output +! +! op_restr - type(psb_dspmat_type), input/output +! The restrictor operator; normally, it is the transpose of the prolongator. +! +! info - integer, output. +! Error code. +! +subroutine amg_d_newmatch_aggregator_inner_mat_asb(ag,parms,a,desc_a,& + & ac,desc_ac, op_prol,op_restr,info) + use psb_base_mod + use amg_base_prec_type +#if defined(SERIAL_MPI) + use amg_d_newmatch_aggregator_mod +#else + use amg_d_newmatch_aggregator_mod, amg_protect_name => amg_d_newmatch_aggregator_inner_mat_asb +#endif + implicit none + class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag + type(amg_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(inout) :: op_prol,op_restr + type(psb_dspmat_type), intent(inout) :: ac + type(psb_desc_type), intent(inout) :: desc_ac + integer(psb_ipk_), intent(out) :: info + ! + type(psb_ctxt_type) :: ictxt + integer(psb_ipk_) :: np, me + type(psb_ld_coo_sparse_mat) :: acoo, bcoo + type(psb_ld_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl, inl + integer(psb_lpk_) :: ntaggr + integer(psb_ipk_) :: err_act, debug_level, debug_unit + character(len=20) :: name='d_newmatch_inner_mat_asb' + character(len=80) :: aname + logical, parameter :: debug=.false., dump_prol_restr=.false. + + + 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) + +#if !defined(SERIAL_MPI) + + if (debug) write(0,*) me,' ',trim(name),' Start:',& + & allocated(ag%ac),allocated(ag%desc_ac), allocated(ag%prol),allocated(ag%restr) + + select case(parms%coarse_mat) + + case(amg_distr_mat_) + ! Do nothing, it has already been done in spmm_bld_ov. + + case(amg_repl_mat_) + ! + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='no repl coarse_mat_ here') + goto 9999 + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid amg_coarse_mat_') + goto 9999 + end select +#endif + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + +end subroutine amg_d_newmatch_aggregator_inner_mat_asb diff --git a/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 index 6a38b7b8..597b39f1 100644 --- a/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 +++ b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 @@ -247,16 +247,16 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& end if if (i == n_sweeps) call tmp_prol%clone(tmp_pg,info) if (do_timings) call psb_tic(idx_spmmbld) -!!$ ! -!!$ ! On entry, prolv(i) is in global numbering, -!!$ ! -!!$ call amg_d_newmatch_spmm_bld_ov(acv(i-1),desc_acv(i-1),ixaggr,nxaggr,parms,& -!!$ & acv(i),desc_acv(i), prolv(i),restrv(1),tmp_prol,info) -!!$ if (psb_errstatus_fatal()) write(0,*)me,trim(name),'Error fatal on exit from bld_ov(i)',info -!!$ if (debug) then -!!$ call psb_barrier(ictxt) -!!$ if (me==0) write(0,*) me,trim(name),' Done spmm_bld:',i -!!$ end if + ! + ! On entry, prolv(i) is in global numbering, + ! + call amg_d_newmatch_spmm_bld_ov(acv(i-1),desc_acv(i-1),ixaggr,nxaggr,parms,& + & acv(i),desc_acv(i), prolv(i),restrv(1),tmp_prol,info) + if (psb_errstatus_fatal()) write(0,*)me,trim(name),'Error fatal on exit from bld_ov(i)',info + if (debug) then + call psb_barrier(ictxt) + if (me==0) write(0,*) me,trim(name),' Done spmm_bld:',i + end if if (do_timings) call psb_toc(idx_spmmbld) ! Keep a copy of prolv(i) in global numbering for the time being, will diff --git a/amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_inner.F90 b/amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_inner.F90 new file mode 100644 index 00000000..982fd00a --- /dev/null +++ b/amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_inner.F90 @@ -0,0 +1,218 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! File: amg_daggrmat_nosmth_bld.F90 +! +! Subroutine: amg_daggrmat_nosmth_bld +! Version: real +! +! This routine builds a coarse-level matrix A_C 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 the piecewise constant interpolation operator corresponding +! the fine-to-coarse level mapping built by amg_aggrmap_bld. +! +! The coarse-level matrix A_C is distributed among the parallel processes or +! replicated on each of them, according to the value of p%parms%coarse_mat +! specified by the user through amg_dprecinit and amg_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 +! +! For 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: +! 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. +! p - type(amg_d_onelev_type), input/output. +! 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. +! parms - type(amg_dml_parms), input +! Parameters controlling the choice of algorithm +! ac - type(psb_dspmat_type), output +! The coarse matrix on output +! +! 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. +! 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 amg_d_newmatch_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,& + & ac,desc_ac,op_prol,op_restr,t_prol,info) + use psb_base_mod + use amg_d_inner_mod +#if defined(SERIAL_MPI) + use amg_d_newmatch_aggregator_mod +#else + use amg_d_newmatch_aggregator_mod, amg_protect_name => amg_d_newmatch_spmm_bld_inner +#endif + implicit none + + ! Arguments + type(psb_d_csr_sparse_mat), intent(inout) :: a_csr + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_ldspmat_type), intent(inout) :: t_prol + type(psb_dspmat_type), intent(inout) :: ac, op_prol, op_restr + type(psb_desc_type), intent(out) :: desc_ac + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + type(psb_ctxt_type) :: ictxt + integer(psb_ipk_) :: np, me, ndx + character(len=40) :: name + type(psb_ld_coo_sparse_mat) :: tmpcoo + type(psb_d_coo_sparse_mat) :: coo_prol, coo_restr + type(psb_d_csr_sparse_mat) :: ac_csr, csr_restr + type(psb_desc_type), target :: tmp_desc + type(psb_ldspmat_type) :: lac + integer(psb_ipk_) :: debug_level, debug_unit, naggr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nrl, nzl, ip, & + & nzt, naggrm1, naggrp1, i, k + integer(psb_lpk_), allocatable :: ia(:),ja(:) + !integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza, nrpsave, ncpsave, nzpsave + logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false. + integer(psb_ipk_), save :: idx_spspmm=-1, idx_prolcnv=-1, idx_proltrans=-1, idx_asb=-1 + + name='amg_newmatch_spmm_bld_inner' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nglob = desc_a%get_global_rows() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + if ((do_timings).and.(idx_spspmm==-1)) & + & idx_spspmm = psb_get_timer_idx("SPMM_BLD: spspmm ") + if ((do_timings).and.(idx_prolcnv==-1)) & + & idx_prolcnv = psb_get_timer_idx("SPMM_BLD: prolcnv ") + if ((do_timings).and.(idx_proltrans==-1)) & + & idx_proltrans = psb_get_timer_idx("SPMM_BLD: proltrans") + if ((do_timings).and.(idx_asb==-1)) & + & idx_asb = psb_get_timer_idx("SPMM_BLD: asb ") + + if (do_timings) call psb_tic(idx_prolcnv) + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + +#if !defined(SERIAL_MPI) + ! + ! Here T_PROL should be arriving with GLOBAL indices on the cols + ! and LOCAL indices on the rows. + ! + if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',& + & op_prol%get_fmt(),op_prol%get_nrows(),op_prol%get_ncols(),op_prol%get_nzeros(),& + & nrow,ntaggr,naggr + + call t_prol%cp_to(tmpcoo) + + call psb_cdall(ictxt,desc_ac,info,nl=naggr) + nzl = tmpcoo%get_nzeros() + if (debug) write(0,*) me,' ',trim(name),' coo_prol: ',& + & tmpcoo%ia(1:min(10,nzl)),' :',tmpcoo%ja(1:min(10,nzl)) + call desc_ac%indxmap%g2lip_ins(tmpcoo%ja(1:nzl),info) + call tmpcoo%set_ncols(desc_ac%get_local_cols()) + call tmpcoo%cp_to_icoo(coo_prol,info) + + call amg_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,& + & coo_prol,desc_ac,coo_restr,info) + + nzl = coo_prol%get_nzeros() + if (debug) write(0,*) me,' ',trim(name),' coo_prol: ',& + & coo_prol%ia(1:min(10,nzl)),' :',coo_prol%ja(1:min(10,nzl)) + + call op_prol%mv_from(coo_prol) + call op_restr%mv_from(coo_restr) + + if (debug) then + write(0,*) me,' ',trim(name),' Checkpoint at exit' + call psb_barrier(ictxt) + write(0,*) me,' ',trim(name),' Checkpoint through' + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x a3') + goto 9999 + end if + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done smooth_aggregate ' +#endif + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine amg_d_newmatch_spmm_bld_inner diff --git a/amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_ov.F90 b/amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_ov.F90 new file mode 100644 index 00000000..8e8fdc53 --- /dev/null +++ b/amgprec/impl/aggregator/amg_d_newmatch_spmm_bld_ov.F90 @@ -0,0 +1,169 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! File: amg_daggrmat_nosmth_bld_ov.F90 +! +! Subroutine: amg_daggrmat_nosmth_bld_ov +! Version: real +! +! This routine builds a coarse-level matrix A_C 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 the piecewise constant interpolation operator corresponding +! the fine-to-coarse level mapping built by amg_aggrmap_bld_ov. +! +! The coarse-level matrix A_C is distributed among the parallel processes or +! replicated on each of them, according to the value of p%parms%coarse_mat +! specified by the user through amg_dprecinit and amg_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 +! +! For 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: +! 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. +! p - type(amg_d_onelev_type), input/output. +! 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. +! parms - type(amg_dml_parms), input +! Parameters controlling the choice of algorithm +! ac - type(psb_dspmat_type), output +! The coarse matrix on output +! +! 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. +! 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 amg_d_newmatch_spmm_bld_ov(a,desc_a,ilaggr,nlaggr,parms,& + & ac,desc_ac,op_prol,op_restr,t_prol,info) + use psb_base_mod + use amg_d_inner_mod +#if defined(SERIAL_MPI) + use amg_d_newmatch_aggregator_mod +#else + use amg_d_newmatch_aggregator_mod, amg_protect_name => amg_d_newmatch_spmm_bld_ov +#endif + implicit none + + ! Arguments + type(psb_dspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_ldspmat_type), intent(inout) :: t_prol + type(psb_dspmat_type), intent(inout) :: ac, op_prol, op_restr + type(psb_desc_type), intent(out) :: desc_ac + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_) :: err_act + + type(psb_ctxt_type) :: ictxt + integer(psb_ipk_) :: np, me + character(len=20) :: name + type(psb_d_csr_sparse_mat) :: acsr + type(psb_ld_coo_sparse_mat) :: coo_prol, coo_restr + integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, & + & naggr, nzt, naggrm1, naggrp1, i, k + integer(psb_ipk_) :: inaggr, nzlp + integer(psb_ipk_) :: debug_level, debug_unit + logical, parameter :: debug=.false., new_version=.true. + + name='amg_newmatch_spmm_bld_ov' + if(psb_get_errstatus().ne.0) return + info=psb_success_ + call psb_erractionsave(err_act) + + + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + +#if !defined(SERIAL_MPI) + + call a%mv_to(acsr) + + call amg_d_newmatch_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,& + & ac,desc_ac,op_prol,op_restr,t_prol,info) + if (psb_errstatus_fatal()) write(0,*)me,trim(name),'Error fatal on exit from bld_inner',info + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err="SPMM_BLD_INNER") + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done spmm_bld ' +#endif + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine amg_d_newmatch_spmm_bld_ov