fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Start work on new matching aggregation

Browse Source
TestFerdous
Salvatore Filippone 3 years ago
parent 9bc8b540b3
commit 2396123b44
7 changed files with 1273 additions and 1 deletions
Show Stats Download Patch File Download Diff File

3
amgprec/Makefile
Unescape Escape View File

@ -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 \

433
amgprec/amg_d_newmatch_aggregator_mod.F90
Unescape Escape View File

@ -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

184
amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_asb.f90
Unescape Escape View File

@ -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

214
amgprec/impl/aggregator/mld_d_bcmatch_aggregator_mat_bld.f90
Unescape Escape View File

@ -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

209
amgprec/impl/aggregator/mld_d_bcmatch_aggregator_tprol.f90
Unescape Escape View File

@ -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

159
amgprec/impl/aggregator/mld_d_bcmatch_map_to_tprol.f90
Unescape Escape View File

@ -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

72
amgprec/impl/aggregator/newmatch_interface.cpp
Unescape Escape View File

@ -0,0 +1,72 @@
#include <string.h>
#include <stdio.h>
#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<irp[i+1]; j++) {
ilaggr[i] = ja[j] + 1;
valaggr[i] = val[j];
}
}
*num_cols = nc;
//w_inp=bcm_VectorData(w);
bcm_CSRMatrixDestroy(Ac);
bcm_CSRMatrixDestroy(P);
return(0);
}
Write Preview
Loading…
Cancel
Save