diff --git a/amgprec/Makefile b/amgprec/Makefile index 0adac3e9..c7438730 100644 --- a/amgprec/Makefile +++ b/amgprec/Makefile @@ -16,7 +16,8 @@ DMODOBJS=amg_d_prec_type.o \ amg_d_dec_aggregator_mod.o amg_d_symdec_aggregator_mod.o \ amg_d_ainv_solver.o amg_d_base_ainv_mod.o \ amg_d_invk_solver.o amg_d_invt_solver.o amg_d_krm_solver.o \ - amg_d_matchboxp_mod.o amg_d_parmatch_aggregator_mod.o + amg_d_matchboxp_mod.o amg_d_parmatch_aggregator_mod.o \ + amd_d_newmatch_aggregator_mod.o SMODOBJS=amg_s_prec_type.o amg_s_ilu_fact_mod.o \ amg_s_inner_mod.o amg_s_ilu_solver.o amg_s_diag_solver.o amg_s_jac_smoother.o amg_s_as_smoother.o \ diff --git a/amgprec/amg_d_newmatch_aggregator_mod.F90 b/amgprec/amg_d_newmatch_aggregator_mod.F90 new file mode 100644 index 00000000..b6c4fc97 --- /dev/null +++ b/amgprec/amg_d_newmatch_aggregator_mod.F90 @@ -0,0 +1,433 @@ +! +! +! The aggregator object hosts the aggregation method for building +! the multilevel hierarchy. This variant is based on the hybrid method +! presented in +! +! +! sm - class(amg_T_base_smoother_type), allocatable +! The current level preconditioner (aka smoother). +! parms - type(amg_RTml_parms) +! The parameters defining the multilevel strategy. +! ac - The local part of the current-level matrix, built by +! coarsening the previous-level matrix. +! desc_ac - type(psb_desc_type). +! The communication descriptor associated to the matrix +! stored in ac. +! base_a - type(psb_Tspmat_type), pointer. +! Pointer (really a pointer!) to the local part of the current +! matrix (so we have a unified treatment of residuals). +! We need this to avoid passing explicitly the current matrix +! to the routine which applies the preconditioner. +! base_desc - type(psb_desc_type), pointer. +! Pointer to the communication descriptor associated to the +! matrix pointed by base_a. +! map - Stores the maps (restriction and prolongation) between the +! vector spaces associated to the index spaces of the previous +! and current levels. +! +! Methods: +! Most methods follow the encapsulation hierarchy: they take whatever action +! is appropriate for the current object, then call the corresponding method for +! the contained object. +! As an example: the descr() method prints out a description of the +! level. It starts by invoking the descr() method of the parms object, +! then calls the descr() method of the smoother object. +! +! descr - Prints a description of the object. +! default - Set default values +! dump - Dump to file object contents +! set - Sets various parameters; when a request is unknown +! it is passed to the smoother object for further processing. +! check - Sanity checks. +! sizeof - Total memory occupation in bytes +! get_nzeros - Number of nonzeros +! +! + +module amg_d_newmatch_aggregator_mod + use amg_d_base_aggregator_mod + use iso_c_binding + + type, bind(c):: bcm_Vector + type(c_ptr) :: data + integer(c_int) :: size + integer(c_int) :: owns_data + end type bcm_Vector + + type, bind(c):: bcm_CSRMatrix + type(c_ptr) :: i + type(c_ptr) :: j + integer(c_int) :: num_rows + integer(c_int) :: num_cols + integer(c_int) :: num_nonzeros + integer(c_int) :: owns_data + type(c_ptr) :: data + end type bcm_CSRMatrix + + type, extends(amg_d_base_aggregator_type) :: amg_d_newmatch_aggregator_type + integer(psb_ipk_) :: matching_alg + integer(psb_ipk_) :: n_sweeps + ! + ! Note: the BootCMatch kernel we invoke overwrites + ! the W argument with its update. Hence, copy it in w_nxt + ! before passing it to the matching + ! + real(psb_dpk_), allocatable :: w(:), w_nxt(:) + type(bcm_Vector) :: w_c_nxt + integer(psb_ipk_) :: max_csize + integer(psb_ipk_) :: max_nlevels + contains + procedure, pass(ag) :: bld_tprol => amg_d_newmatch_aggregator_build_tprol + procedure, pass(ag) :: cseti => d_newmatch_aggr_cseti + procedure, pass(ag) :: default => d_newmatch_aggr_set_default + procedure, pass(ag) :: mat_asb => amg_d_newmatch_aggregator_mat_asb + procedure, pass(ag) :: mat_bld => amg_d_newmatch_aggregator_mat_bld + procedure, pass(ag) :: update_next => d_newmatch_aggregator_update_next + procedure, pass(ag) :: bld_wnxt => d_newmatch_bld_wnxt + procedure, pass(ag) :: bld_default_w => d_bld_default_w + procedure, pass(ag) :: set_c_default_w => d_set_default_bcm_w + procedure, pass(ag) :: descr => d_newmatch_aggregator_descr + procedure, pass(ag) :: clone => d_newmatch_aggregator_clone + procedure, pass(ag) :: free => d_newmatch_aggregator_free + procedure, nopass :: fmt => d_newmatch_aggregator_fmt + end type amg_d_newmatch_aggregator_type + + + interface + subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,t_prol,info) + import :: amg_d_newmatch_aggregator_type, psb_desc_type, & + & psb_dspmat_type, psb_ldspmat_type, psb_dpk_, & + & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms, amg_daggr_data + implicit none + class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag + type(amg_dml_parms), intent(inout) :: parms + type(amg_daggr_data), intent(in) :: ag_data + type(psb_dspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_ldspmat_type), intent(out) :: t_prol + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_newmatch_aggregator_build_tprol + end interface + + interface + subroutine amg_d_newmatch_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& + & 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_, psb_epk_, amg_dml_parms + 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(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol,ac,op_restr + type(psb_ldspmat_type), intent(inout) :: t_prol + type(psb_desc_type), intent(inout) :: desc_ac + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_newmatch_aggregator_mat_bld + end interface + + interface + subroutine amg_d_newmatch_aggregator_mat_asb(ag,parms,a,desc_a,& + & ac,desc_ac, op_prol,op_restr,info) + import :: amg_d_newmatch_aggregator_type, psb_desc_type, & + & psb_dspmat_type, psb_ldspmat_type, psb_dpk_, & + & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms + 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(inout) :: desc_a + type(psb_dspmat_type), intent(inout) :: op_prol,ac,op_restr + type(psb_desc_type), intent(inout) :: desc_ac + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_newmatch_aggregator_mat_asb + end interface + + + interface + subroutine amg_d_newmatch_map_to_tprol(desc_a,ilaggr,nlaggr,valaggr, op_prol,info) + import :: amg_d_newmatch_aggregator_type, psb_desc_type, & + & psb_dspmat_type, psb_ldspmat_type, psb_dpk_, & + & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:) + real(psb_dpk_), allocatable, intent(inout) :: valaggr(:) + type(psb_ldspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_newmatch_map_to_tprol + end interface + + interface + subroutine amg_daggrmat_unsmth_spmm_asb(a,desc_a,ilaggr,nlaggr,parms,& + & ac,op_prol,op_restr,info) + import :: amg_d_newmatch_aggregator_type, psb_desc_type, & + & psb_dspmat_type, psb_ldspmat_type, psb_dpk_, & + & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms + implicit none + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(amg_dml_parms), intent(inout) :: parms + type(psb_ldspmat_type), intent(inout) :: op_prol + type(psb_ldspmat_type), intent(out) :: ac,op_restr + integer(psb_ipk_), intent(out) :: info + end subroutine amg_daggrmat_unsmth_spmm_asb + end interface + +!!$ interface +!!$ subroutine amg_d_newmatch_unsmth_spmm_bld(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_, psb_epk_, amg_dml_parms +!!$ +!!$ implicit none +!!$ +!!$ ! Arguments +!!$ type(psb_dspmat_type), intent(in) :: a +!!$ type(psb_desc_type), intent(in) :: 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(inout) :: desc_ac +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine amg_d_newmatch_unsmth_spmm_bld +!!$ end interface + + + private :: is_legal_malg, is_legal_csize, is_legal_nsweeps, is_legal_nlevels + + +contains + + subroutine d_bld_default_w(ag,nr) + use psb_realloc_mod + implicit none + class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag + integer(psb_ipk_), intent(in) :: nr + integer(psb_ipk_) :: info + call psb_realloc(nr,ag%w,info) + if (info /= psb_success_) return + ag%w = done + call ag%set_c_default_w() + end subroutine d_bld_default_w + + subroutine d_set_default_bcm_w(ag) + use psb_realloc_mod + use iso_c_binding + implicit none + class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag + integer(psb_ipk_) :: info + + call psb_safe_ab_cpy(ag%w,ag%w_nxt,info) + ag%w_c_nxt%size = psb_size(ag%w_nxt) + ag%w_c_nxt%owns_data = 0 + if (ag%w_c_nxt%size > 0) call set_cloc(ag%w_nxt, ag%w_c_nxt) + + end subroutine d_set_default_bcm_w + + subroutine set_cloc(vect,w_c_nxt) + use iso_c_binding + real(psb_dpk_), target :: vect(:) + type(bcm_Vector) :: w_c_nxt + + w_c_nxt%data = c_loc(vect) + end subroutine set_cloc + + + subroutine d_newmatch_bld_wnxt(ag,ilaggr,valaggr,nx) + use psb_realloc_mod + implicit none + class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag + integer(psb_lpk_), intent(in) :: ilaggr(:) + real(psb_dpk_), intent(in) :: valaggr(:) + integer(psb_ipk_), intent(in) :: nx + + integer(psb_ipk_) :: info,i,j + + ! The vector was already fixed in the call to Newmatch. + call psb_realloc(nx,ag%w_nxt,info) + + end subroutine d_newmatch_bld_wnxt + + function d_newmatch_aggregator_fmt() result(val) + implicit none + character(len=32) :: val + + val = "BootCMatch aggregation" + end function d_newmatch_aggregator_fmt + + subroutine d_newmatch_aggregator_descr(ag,parms,iout,info) + implicit none + class(amg_d_newmatch_aggregator_type), intent(in) :: ag + type(amg_dml_parms), intent(in) :: parms + integer(psb_ipk_), intent(in) :: iout + integer(psb_ipk_), intent(out) :: info + + write(iout,*) 'BootCMatch Aggregator' + write(iout,*) ' Number of BootCMatch sweeps: ',ag%n_sweeps + write(iout,*) ' Matching algorithm : ',ag%matching_alg + write(iout,*) ' 0: Preis 1: MC64 2: SPRAL ' + write(iout,*) 'Aggregator object type: ',ag%fmt() + call parms%mldescr(iout,info) + + return + end subroutine d_newmatch_aggregator_descr + + function is_legal_malg(alg) result(val) + logical :: val + integer(psb_ipk_) :: alg + + val = ((0<=alg).and.(alg<=2)) + end function is_legal_malg + + function is_legal_csize(csize) result(val) + logical :: val + integer(psb_ipk_) :: csize + + val = ((-1==csize).or.(csize >0)) + end function is_legal_csize + + function is_legal_nsweeps(nsw) result(val) + logical :: val + integer(psb_ipk_) :: nsw + + val = (1<=nsw) + end function is_legal_nsweeps + + function is_legal_nlevels(nlv) result(val) + logical :: val + integer(psb_ipk_) :: nlv + + val = (1<=nlv) + end function is_legal_nlevels + + subroutine d_newmatch_aggregator_update_next(ag,agnext,info) + use psb_realloc_mod + implicit none + class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag + class(amg_d_base_aggregator_type), target, intent(inout) :: agnext + integer(psb_ipk_), intent(out) :: info + + ! + ! + select type(agnext) + class is (amg_d_newmatch_aggregator_type) + if (.not.is_legal_malg(agnext%matching_alg)) & + & agnext%matching_alg = ag%matching_alg + if (.not.is_legal_nsweeps(agnext%n_sweeps))& + & agnext%n_sweeps = ag%n_sweeps + if (.not.is_legal_csize(agnext%max_csize))& + & agnext%max_csize = ag%max_csize + if (.not.is_legal_nlevels(agnext%max_nlevels))& + & agnext%max_nlevels = ag%max_nlevels + ! Is this going to generate shallow copies/memory leaks/double frees? + ! To be investigated further. + call psb_safe_ab_cpy(ag%w_nxt,agnext%w,info) + call agnext%set_c_default_w() + class default + ! What should we do here? + end select + info = 0 + end subroutine d_newmatch_aggregator_update_next + + subroutine d_newmatch_aggr_cseti(ag,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_d_newmatch_aggregator_type), intent(inout) :: ag + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act, iwhat + character(len=20) :: name='d_newmatch_aggr_cseti' + info = psb_success_ + + ! For now we ignore IDX + + select case(what) + case('BCM_MATCH_ALG') + ag%matching_alg=val + case('BCM_SWEEPS') + ag%n_sweeps=val + case('BCM_MAX_CSIZE') + ag%max_csize=val + case('BCM_MAX_NLEVELS') + ag%max_nlevels=val + case('BCM_W_SIZE') + call ag%bld_default_w(val) + case default + ! Do nothing + end select + return + end subroutine d_newmatch_aggr_cseti + + subroutine d_newmatch_aggr_set_default(ag) + + Implicit None + + ! Arguments + class(amg_d_newmatch_aggregator_type), intent(inout) :: ag + character(len=20) :: name='d_newmatch_aggr_set_default' + ag%matching_alg = 0 + ag%n_sweeps = 1 + ag%max_nlevels = 36 + ag%max_csize = -1 + ! + ! Apparently BootCMatch works better + ! by keeping all entries + ! + ag%do_clean_zeros = .false. + + return + + end subroutine d_newmatch_aggr_set_default + + subroutine d_newmatch_aggregator_free(ag,info) + use iso_c_binding + implicit none + class(amg_d_newmatch_aggregator_type), intent(inout) :: ag + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(ag%w)) deallocate(ag%w,stat=info) + if (info /= 0) return + if (allocated(ag%w_nxt)) deallocate(ag%w_nxt,stat=info) + if (info /= 0) return + ag%w_c_nxt%size = 0 + ag%w_c_nxt%data = c_null_ptr + ag%w_c_nxt%owns_data = 0 + end subroutine d_newmatch_aggregator_free + + subroutine d_newmatch_aggregator_clone(ag,agnext,info) + implicit none + class(amg_d_newmatch_aggregator_type), intent(inout) :: ag + class(amg_d_base_aggregator_type), allocatable, intent(inout) :: agnext + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(agnext)) then + call agnext%free(info) + if (info == 0) deallocate(agnext,stat=info) + end if + if (info /= 0) return + allocate(agnext,source=ag,stat=info) + select type(agnext) + class is (amg_d_newmatch_aggregator_type) + call agnext%set_c_default_w() + class default + ! Should never ever get here + info = -1 + end select + end subroutine d_newmatch_aggregator_clone + +end module amg_d_newmatch_aggregator_mod diff --git a/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_asb.f90 b/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_asb.f90 new file mode 100644 index 00000000..92609a14 --- /dev/null +++ b/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_asb.f90 @@ -0,0 +1,184 @@ +! +! +! MLD2P4 version 2.2 +! MultiLevel Domain Decomposition Parallel Preconditioners Package +! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the 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_mat_asb.f90 +! +! Subroutine: mld_d_bcmatch_aggregator_mat_asb +! Version: real +! +! +! From a given AC to final format, generating DESC_AC +! +! Arguments: +! ag - type(mld_d_bcmatch_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), 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 mld_d_bcmatch_aggregator_mat_asb(ag,parms,a,desc_a,& + & ac,desc_ac, op_prol,op_restr,info) + use psb_base_mod + use mld_base_prec_type + use mld_d_bcmatch_aggregator_mod, mld_protect_name => mld_d_bcmatch_aggregator_mat_asb + 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(inout) :: desc_a + type(psb_dspmat_type), intent(inout) :: op_prol, ac,op_restr + type(psb_desc_type), intent(inout) :: desc_ac + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: ictxt, np, me + type(psb_ld_coo_sparse_mat) :: tmpcoo + type(psb_ldspmat_type) :: tmp_ac + integer(psb_ipk_) :: i_nr, i_nc, i_nl, nzl + integer(psb_lpk_) :: ntaggr + integer(psb_ipk_) :: err_act, debug_level, debug_unit + character(len=20) :: name='d_bcmatch_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) + + select case(parms%coarse_mat) + + case(mld_distr_mat_) + + call ac%cscnv(info,type='csr') + call op_prol%cscnv(info,type='csr') + call op_restr%cscnv(info,type='csr') + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ac ' + + case(mld_repl_mat_) + ! + ! We are assuming here that an d matrix + ! can hold all entries + ! + if (desc_ac%get_global_rows() < huge(1_psb_ipk_) ) then + ntaggr = desc_ac%get_global_rows() + i_nr = ntaggr + else + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end if + + call op_prol%mv_to(tmpcoo) + nzl = tmpcoo%get_nzeros() + call psb_loc_to_glob(tmpcoo%ja(1:nzl),desc_ac,info,'I') + call op_prol%mv_from(tmpcoo) + + call op_restr%mv_to(tmpcoo) + nzl = tmpcoo%get_nzeros() + call psb_loc_to_glob(tmpcoo%ia(1:nzl),desc_ac,info,'I') + call op_restr%mv_from(tmpcoo) + + call op_prol%set_ncols(i_nr) + call op_restr%set_nrows(i_nr) + + call psb_gather(tmp_ac,ac,desc_ac,info,root=-ione,& + & dupl=psb_dupl_add_,keeploc=.false.) + call tmp_ac%mv_to(tmpcoo) + call ac%mv_from(tmpcoo) + + call psb_cdall(ictxt,desc_ac,info,mg=ntaggr,repl=.true.) + if (info == psb_success_) call psb_cdasb(desc_ac,info) + ! + ! Now that we have the descriptors and the restrictor, we should + ! update the W. But we don't, because REPL is only valid + ! at the coarsest level, so no need to carry over. + ! + + if (info /= psb_success_) goto 9999 + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end select + + 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/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_bld.f90 b/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_bld.f90 new file mode 100644 index 00000000..6ea944c1 --- /dev/null +++ b/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_bld.f90 @@ -0,0 +1,214 @@ +! +! +! MLD2P4 Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, 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_bld.f90 +! +! Subroutine: mld_d_base_aggregator_mat_bld +! 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_bld. +! +! 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_bld(ag,parms,a,desc_a,ilaggr,nlaggr,& + & ac,desc_ac,op_prol,op_restr,t_prol,info) + use psb_base_mod + use mld_d_inner_mod + use mld_d_prec_type + use mld_d_bcmatch_aggregator_mod, mld_protect_name => mld_d_bcmatch_aggregator_mat_bld + 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(inout) :: desc_a + integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(psb_ldspmat_type), intent(inout) :: t_prol + type(psb_dspmat_type), intent(inout) :: op_prol,ac,op_restr + type(psb_desc_type), intent(inout) :: desc_ac + integer(psb_ipk_), intent(out) :: info + + ! Local variables + character(len=20) :: name + integer(psb_mpk_) :: ictxt, np, me + type(psb_ld_coo_sparse_mat) :: acoo, bcoo + type(psb_ld_csr_sparse_mat) :: acsr1 + integer(psb_lpk_) :: nzl,ntaggr + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_d_bcmatch_aggregator_mat_bld' + 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) + + ! + ! 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 + ! + select case (parms%aggr_prol) + case (mld_no_smooth_) + +!!$ call mld_d_bcmatch_unsmth_spmm_bld(a,desc_a,ilaggr,nlaggr,& +!!$ & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + call mld_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr, & + & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + + + case(mld_smooth_prol_) + + call mld_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, & + & parms,ac,desc_ac,op_prol,op_restr,t_prol,info) + +!!$ case(mld_biz_prol_) +!!$ +!!$ call mld_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, & +!!$ & parms,ac,desc_ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, & + & parms,ac,desc_ac,op_prol,op_restr,t_prol,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 bld') + 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_bld diff --git a/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_tprol.f90 b/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_tprol.f90 new file mode 100644 index 00000000..941dd05e --- /dev/null +++ b/amgprec/impl/aggregator/mld_d_bcmatch_aggregator_tprol.f90 @@ -0,0 +1,209 @@ +! +! +! MLD2P4 Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, 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 +! +! + +subroutine mld_d_bcmatch_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod + use mld_d_prec_type + use mld_d_bcmatch_aggregator_mod, mld_protect_name => mld_d_bcmatch_aggregator_build_tprol + use mld_d_inner_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(mld_daggr_data), intent(in) :: ag_data + type(psb_dspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a + integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_ldspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + + + ! Local variables + real(psb_dpk_), allocatable:: valaggr(:) + 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_mpk_) :: 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(:), c_ilaggr(:) + 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 + import + 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 + import + 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_cycle,'Multilevel cycle',& + & mld_mult_ml_,is_legal_ml_cycle) + call mld_check_def(parms%par_aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_par_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) + if (ag%do_clean_zeros) call acsr%clean_zeros(info) + nr = a%get_nrows() + if (psb_size(ag%w) < nr) call ag%bld_default_w(nr) + + !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,c_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 + if (ag%max_csize > 0) then + max_csize = ag%max_csize + else + max_csize = (ag_data%min_coarse_size + np -1)/np + end if + if (ag%max_nlevels > 0) then + max_nlevels = ag%max_nlevels + else + max_nlevels = ag_data%max_levs + end if + + info = mld_bootCMatch_if(C,match_algorithm,n_sweeps,max_nlevels,max_csize,& + & ag%w_c_nxt, isz, c_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 + ilaggr(1:nr) = c_ilaggr(1:nr) +!!$ write(0,*) 'On output from BootCMatch',nr,num_pcols,size(ilaggr),maxval(ilaggr),& +!!$ & minval(ilaggr),minval(ilaggr(1:nr)),a%get_nrows(),a%get_ncols() + ! Prepare vector W for next level, just in case + call ag%bld_wnxt(ilaggr(1:nr),valaggr(1:nr),num_pcols) + + 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_d_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/amgprec/impl/aggregator/mld_d_bcmatch_map_to_tprol.f90 b/amgprec/impl/aggregator/mld_d_bcmatch_map_to_tprol.f90 new file mode 100644 index 00000000..2ab92c4f --- /dev/null +++ b/amgprec/impl/aggregator/mld_d_bcmatch_map_to_tprol.f90 @@ -0,0 +1,159 @@ +! +! +! MLD2P4 Extensions +! +! (C) Copyright 2019 +! +! Salvatore Filippone Cranfield University +! Pasqua D'Ambra IAC-CNR, Naples, 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 + use mld_d_bcmatch_aggregator_mod, mld_protect_name => mld_d_bcmatch_map_to_tprol + + implicit none + + ! Arguments + type(psb_desc_type), intent(in) :: desc_a + integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:) + real(psb_dpk_), allocatable, intent(inout) :: valaggr(:) + type(psb_ldspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_lpk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr + type(psb_ld_coo_sparse_mat) :: tmpcoo + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_lpk_) :: 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) + j = 0 + do i=1,ncol + if (valaggr(i) /= dzero) then + j = j + 1 + tmpcoo%val(j) = valaggr(i) + tmpcoo%ia(j) = i + tmpcoo%ja(j) = ilaggr(i) + end if + end do + call tmpcoo%set_nzeros(j) + 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/amgprec/impl/aggregator/newmatch_interface.cpp b/amgprec/impl/aggregator/newmatch_interface.cpp new file mode 100644 index 00000000..d02dd398 --- /dev/null +++ b/amgprec/impl/aggregator/newmatch_interface.cpp @@ -0,0 +1,72 @@ + +#include +#include + +#include "bcm.h" + +bcm_CSRMatrix bootCMatch(bcm_CSRMatrix *C, int *match_algorithm, int *n_sweeps, int *max_nlevels, int *max_csize, bcm_Vector *w); +bcm_CSRMatrix bootCMatch(bcm_CSRMatrix *C, int *match_algorithm, int *n_sweeps, int *max_nlevels, int *max_csize, bcm_Vector *w){ + bcm_Vector *w_temp; + int info; + //double *w_inp; + //w_inp=bcm_VectorData(w); + + bcm_CSRMatrix *P; + bcm_CSRMatrix *Ac; + int ftcoarse=1; + int cr_it=0, cr_relax_type=0; + double cr_relax_weight=0.0; + // Here I am building Ac but I won't use it. + Ac=bcm_CSRMatchingAgg(C, &w, &P, *match_algorithm, *n_sweeps,*max_csize, *max_nlevels, &ftcoarse, + cr_it, cr_relax_type, cr_relax_weight); + //w_inp=bcm_VectorData(w); + bcm_CSRMatrixDestroy(Ac); + return *P; +} + +int mld_bootCMatch_if(bcm_CSRMatrix *C, int match_algorithm, int n_sweeps, + int max_nlevels, int max_csize, bcm_Vector *w, + int isz, int ilaggr[], double valaggr[], int *num_cols){ + bcm_Vector *w_temp; + int info; + //double *w_inp; + //w_inp=bcm_VectorData(w); + + bcm_CSRMatrix *P; + bcm_CSRMatrix *Ac; + int *irp, *ja, nr, nz, nc,i,j; + double *val; + int ftcoarse=1; + int cr_it=0, cr_relax_type=0; + double cr_relax_weight=0.0; + + // Sanity checks + nr = bcm_CSRMatrixNumRows(C); + nc = bcm_VectorSize(w); +// fprintf(stderr,"Sanity check: %d %d \n",nr,nc); + + // Here I am building Ac but I won't use it. + Ac=bcm_CSRMatchingAgg(C, &w, &P, match_algorithm, n_sweeps, max_csize, max_nlevels, &ftcoarse, + cr_it, cr_relax_type, cr_relax_weight); + irp = bcm_CSRMatrixI(P); + ja = bcm_CSRMatrixJ(P); + val = bcm_CSRMatrixData(P); + nr = bcm_CSRMatrixNumRows(P); + nc = bcm_CSRMatrixNumCols(P); + nz = bcm_CSRMatrixNumNonzeros(P); + + if (isz < nr) return(-1); + if (nz != nr) return(-2); + /* loop here only makes sense when nr==nz */ + for (i=0; i< nr; i++) { + for (j=irp[i]; j