From 24988bcfc5ac1385d852741352f561c24d273731 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sat, 5 May 2018 14:35:47 +0100 Subject: [PATCH] Moved BootCMatch interface to tests; to be reviewed later. --- mlprec/Makefile | 3 +- .../mld_d_bcmatch_aggregator_mat_asb.f90 | 225 +++++++++++ .../Bcmatch}/mld_d_bcmatch_aggregator_mod.F90 | 0 .../mld_d_bcmatch_aggregator_tprol.f90 | 355 ++++++++++++++++++ tests/Bcmatch/mld_d_bcmatch_map_bld.f90 | 0 tests/Bcmatch/mld_d_bcmatch_map_to_tprol.f90 | 158 ++++++++ tests/pdegen/runs/mld_pde3d.inp | 2 +- 7 files changed, 741 insertions(+), 2 deletions(-) create mode 100644 tests/Bcmatch/mld_d_bcmatch_aggregator_mat_asb.f90 rename {mlprec => tests/Bcmatch}/mld_d_bcmatch_aggregator_mod.F90 (100%) create mode 100644 tests/Bcmatch/mld_d_bcmatch_aggregator_tprol.f90 create mode 100644 tests/Bcmatch/mld_d_bcmatch_map_bld.f90 create mode 100644 tests/Bcmatch/mld_d_bcmatch_map_to_tprol.f90 diff --git a/mlprec/Makefile b/mlprec/Makefile index d89347b7..5de506ee 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -13,7 +13,8 @@ DMODOBJS=mld_d_prec_type.o mld_d_ilu_fact_mod.o \ mld_d_base_solver_mod.o mld_d_base_smoother_mod.o mld_d_onelev_mod.o \ mld_d_gs_solver.o mld_d_mumps_solver.o \ mld_d_base_aggregator_mod.o mld_d_hybrid_aggregator_mod.o \ - mld_d_symdec_aggregator_mod.o mld_d_bcmatch_aggregator_mod.o + mld_d_symdec_aggregator_mod.o +#mld_d_bcmatch_aggregator_mod.o SMODOBJS=mld_s_prec_type.o mld_s_ilu_fact_mod.o \ mld_s_inner_mod.o mld_s_ilu_solver.o mld_s_diag_solver.o mld_s_jac_smoother.o mld_s_as_smoother.o \ diff --git a/tests/Bcmatch/mld_d_bcmatch_aggregator_mat_asb.f90 b/tests/Bcmatch/mld_d_bcmatch_aggregator_mat_asb.f90 new file mode 100644 index 00000000..750c068b --- /dev/null +++ b/tests/Bcmatch/mld_d_bcmatch_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_bcmatch_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_unsmth_spmm_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_bcmatch_aggregator_mat_asb diff --git a/mlprec/mld_d_bcmatch_aggregator_mod.F90 b/tests/Bcmatch/mld_d_bcmatch_aggregator_mod.F90 similarity index 100% rename from mlprec/mld_d_bcmatch_aggregator_mod.F90 rename to tests/Bcmatch/mld_d_bcmatch_aggregator_mod.F90 diff --git a/tests/Bcmatch/mld_d_bcmatch_aggregator_tprol.f90 b/tests/Bcmatch/mld_d_bcmatch_aggregator_tprol.f90 new file mode 100644 index 00000000..d6b511fe --- /dev/null +++ b/tests/Bcmatch/mld_d_bcmatch_aggregator_tprol.f90 @@ -0,0 +1,355 @@ +! +! +! 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_bcmatch_aggregator_tprol.f90 +! +! Subroutine: mld_d_bcmatch_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. +! + +module bcm_CSRMatrix_mod +use psb_base_mod +use psb_util_mod +use iso_c_binding +use bcm_csr_type_mod +implicit none + +contains + subroutine MLD_to_CSR(a,csr_ia, csr_ja, csr_val, C, info) + type(psb_dspmat_type), intent(in) :: a + type(bcm_CSRMatrix), intent(out) :: C + real(c_double), allocatable, target, intent(out) :: csr_val(:) + integer(c_int), allocatable, target, intent(out) :: csr_ia(:), csr_ja(:) + + !Local variable + character(len=20) :: name + integer(psb_ipk_) :: info + real(psb_dpk_), allocatable :: coo_val(:) + integer(psb_ipk_), allocatable :: coo_ia(:), coo_ja(:) + integer(psb_ipk_) :: x, nz, num_rows, num_cols , i , j, iad, k , k0 + type(psb_d_csr_sparse_mat) :: acsr + + name="MLD_to_CSR" + num_rows= a%get_nrows() + num_cols= a%get_ncols() + nz= a%get_nzeros() + call a%csgetrow(1,num_rows,nz,coo_ia,coo_ja,coo_val,info) + !allocate(csr_ia(0:nz-1), csr_ja(0:nz-1), csr_val(0:nz-1), STAT=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/nz,izero,izero,izero,izero/),& + & a_err='integer') + return + end if + + call a%cp_to(acsr) + allocate(csr_ia(0:nz-1), csr_ja(0:nz-1), csr_val(0:nz-1), STAT=info) + + + csr_ia(0:min(nz,size(acsr%irp,1))-1)=acsr%irp(1:min(nz,size(acsr%irp,1)))-1 + csr_ja(0:nz-1)=acsr%ja(1:nz)-1 + csr_val(0:nz-1)=acsr%val(1:nz) + + call acsr%free() + + C%num_rows=num_rows + C%num_cols=num_cols + C%num_nonzeros=nz + C%owns_data=0 + C%i=c_loc(csr_ia) + C%j=c_loc(csr_ja) + C%data=c_loc(csr_val) + end subroutine MLD_to_CSR + + subroutine bcm_to_op_prol(P, ilaggr, valaggr, info) + type(bcm_CSRMatrix), intent(in) :: P + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:) + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_), allocatable, intent(out) :: valaggr(:) + + ! Local variables + integer(psb_ipk_), pointer :: point_ia(:), point_ja(:) + real(psb_dpk_), pointer :: point_val(:) + integer(psb_ipk_) :: i, j, k + integer(psb_ipk_) :: n, num_rows, num_cols, num_nz + character(len=20) :: name + integer(c_int), allocatable, target :: coo_ia(:),coo_ja(:) + real(c_double), allocatable, target :: coo_val(:) + + name="bcm_to_op_prol" + num_nz = P%num_nonzeros + num_rows = P%num_rows + num_cols = P%num_cols + call c_f_pointer(P%i,point_ia,(/num_rows+1/)) + call c_f_pointer(P%j,point_ja,(/num_nz/)) + call c_f_pointer(P%data,point_val,(/num_nz/)) + + !These are I, J, VAL. + !These are I, J, VAL. + if (allocated(coo_ia)) deallocate(coo_ia) + if (allocated(coo_ja)) deallocate(coo_ja) + if (allocated(coo_val)) deallocate(coo_val) + allocate(coo_ia(num_nz),coo_ja(num_nz),coo_val(num_nz), STAT=info) +!!$ write(0,*) num_rows,num_cols,num_nz,info,size(coo_val) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/num_nz,izero,izero,izero,izero/),& + & a_err='integer') + return + end if + + n = 1 + !coo_ia=-123 + !coo_ja=-123 + coo_val(1:num_nz)=point_val(1:num_nz) + do i=1, num_rows + do j=point_ia(i)+1, point_ia(i+1) + coo_ia(n)=i + coo_ja(n)=point_ja(j) + 1 + n = n + 1 + enddo + enddo + if (allocated(ilaggr)) deallocate(ilaggr) + if (allocated(valaggr)) deallocate(valaggr) + allocate(ilaggr(num_rows),valaggr(num_rows), STAT=info) + ilaggr=0 + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/num_rows,izero,izero,izero,izero/),& + & a_err='integer') + return + end if + + do k=1,num_nz + i=coo_ia(k) + j=coo_ja(k) + ilaggr(i)=j + valaggr(i)=coo_val(i) + enddo + if (allocated(coo_ia)) deallocate(coo_ia) + if (allocated(coo_ja)) deallocate(coo_ja) + if (allocated(coo_val)) deallocate(coo_val) + nullify(point_ia) + nullify(point_ja) + nullify(point_val) + end subroutine bcm_to_op_prol +end module bcm_CSRMatrix_mod + + +subroutine mld_d_bcmatch_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_d_bcmatch_aggregator_mod, mld_protect_name => mld_d_bcmatch_aggregator_build_tprol + use mld_d_inner_mod + !use bcm_CSRMatrix_mod + use bcm_csr_type_mod + use iso_c_binding + implicit none + class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + real(psb_dpk_), allocatable:: valaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + + + ! Local variables + type(psb_dspmat_type) :: a_tmp + type(bcm_CSRMatrix) :: C, P + integer(c_int) :: match_algorithm, n_sweeps, max_csize, max_nlevels + character(len=20) :: name, ch_err + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act, ierr + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: i, j, k, nr, nc, isz, num_pcols + type(psb_d_csr_sparse_mat), target :: acsr + integer(psb_ipk_), allocatable, target :: csr_ia(:), csr_ja(:) + integer(psb_ipk_), allocatable :: aux(:) + real(psb_dpk_), allocatable, target:: csr_val(:) + interface + function bootCMatch(C,match_alg,n_sweeps,max_nlevels,max_csize,w)bind(c,name='bootCMatch') result(P) + use iso_c_binding + use bcm_CSRMatrix_mod + implicit none + type(bcm_CSRMatrix) :: C, P + type(bcm_Vector) :: w + integer(c_int) :: match_alg + integer(c_int) :: n_sweeps + integer(c_int) :: max_nlevels + integer(c_int) :: max_csize + end function bootCMatch + end interface + + interface + function mld_bootCMatch_if(C,match_alg,n_sweeps,max_nlevels,max_csize,& + & w,isz,ilaggr,valaggr, num_cols) & + & bind(c,name='mld_bootCMatch_if') result(iret) + use iso_c_binding + use bcm_CSRMatrix_mod + implicit none + type(bcm_CSRMatrix) :: C, P + type(bcm_Vector) :: w + integer(c_int), value :: match_alg + integer(c_int), value :: n_sweeps + integer(c_int), value :: max_nlevels + integer(c_int), value :: max_csize + integer(c_int), value :: isz + integer(c_int) :: num_cols + integer(c_int) :: ilaggr(*) + real(c_double) :: valaggr(*) + integer(c_int) :: iret + end function mld_bootCMatch_if + end interface + + name='mld_d_bcmatch_aggregator_tprol' + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + 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_ + + + 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 a%csclip(b=a_tmp, info=info, jmax=a%get_nrows(), imax=a%get_nrows()) + + call a_tmp%mv_to(acsr) + + !write(*,*) 'Build_tprol:',acsr%get_nrows(),acsr%get_ncols() + C%num_rows=acsr%get_nrows() + C%num_cols=acsr%get_ncols() + C%num_nonzeros=acsr%get_nzeros() + C%owns_data=0 + acsr%irp = acsr%irp - 1 + acsr%ja = acsr%ja - 1 + C%i=c_loc(acsr%irp) + C%j=c_loc(acsr%ja) + C%data=c_loc(acsr%val) + + isz = a%get_ncols() + call psb_realloc(isz,ilaggr,info) + if (info == psb_success_) call psb_realloc(isz,valaggr,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 + + match_algorithm = ag%matching_alg + n_sweeps = ag%n_sweeps + max_csize = ag%max_csize + max_nlevels = ag%max_nlevels + + info = mld_bootCMatch_if(C,match_algorithm,n_sweeps,max_nlevels,max_csize,& + & ag%w_par, isz, ilaggr, valaggr, num_pcols) + if (info /= psb_success_) then +!!$ write(0,*) 'On return from bootCMatch_if:',info + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_bootCMatch_if') + goto 9999 + end if + + call psb_realloc(np,nlaggr,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 + call acsr%free() + + nlaggr(:)=0 + nlaggr(me+1) = num_pcols + call psb_sum(ictxt,nlaggr(1:np)) + + + call mld_bcmatch_map_to_tprol(desc_a,ilaggr,nlaggr,valaggr,op_prol,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_bcmatch_map_to_tprol') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_d_bcmatch_aggregator_build_tprol diff --git a/tests/Bcmatch/mld_d_bcmatch_map_bld.f90 b/tests/Bcmatch/mld_d_bcmatch_map_bld.f90 new file mode 100644 index 00000000..e69de29b diff --git a/tests/Bcmatch/mld_d_bcmatch_map_to_tprol.f90 b/tests/Bcmatch/mld_d_bcmatch_map_to_tprol.f90 new file mode 100644 index 00000000..960f9cb9 --- /dev/null +++ b/tests/Bcmatch/mld_d_bcmatch_map_to_tprol.f90 @@ -0,0 +1,158 @@ +! +! +! 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_bcmatch_map_to_tprol.f90 +! +! Subroutine: mld_d_bcmatch_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_bcmatch_map_to_tprol(desc_a,ilaggr,nlaggr,valaggr, op_prol,info) + + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_bcmatch_map_to_tprol + + implicit none + + ! Arguments + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:) + real(psb_dpk_), allocatable, intent(inout) :: valaggr(:) + 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_d_bcmatch_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 psb_halo(valaggr,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) = valaggr(i) + 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_bcmatch_map_to_tprol diff --git a/tests/pdegen/runs/mld_pde3d.inp b/tests/pdegen/runs/mld_pde3d.inp index 316edc63..7d1bcbb7 100644 --- a/tests/pdegen/runs/mld_pde3d.inp +++ b/tests/pdegen/runs/mld_pde3d.inp @@ -1,7 +1,7 @@ %%%%%%%%%%% General arguments % Lines starting with % are ignored. CSR ! Storage format CSR COO JAD 0080 ! IDIM; domain size. Linear system size is IDIM**3 -CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES +BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC 00500 ! ITMAX 1 ! ITRACE