Merge branch 'maint-1.0' into development
Salvatore Filippone
3 years ago
commit
b387308455
Binary file not shown.
Binary file not shown.
|
Before Width: | Height: | Size: 905 B After Width: | Height: | Size: 892 B |
Binary file not shown.
|
Before Width: | Height: | Size: 2.1 KiB After Width: | Height: | Size: 2.0 KiB |
File diff suppressed because it is too large
Load Diff
@ -1,59 +0,0 @@
|
||||
MLDDIR=../..
|
||||
MLDINCDIR=$(MLDDIR)/include
|
||||
include $(MLDINCDIR)/Make.inc.amg4psblas
|
||||
MLDMODDIR=$(MLDDIR)/modules
|
||||
MLDLIBDIR=$(MLDDIR)/lib
|
||||
MLD_LIBS=-L$(MLDLIBDIR) -lpsb_krylov -lmld_prec -lpsb_prec
|
||||
FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDMODDIR) $(FMFLAG)$(MLDINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
|
||||
|
||||
HSL_DIR=/opt/hsl/2.3.1/gnu/6.4.0
|
||||
HSL_INCDIR=$(HSL_DIR)/include
|
||||
HSL_LIBDIR=$(HSL_DIR)/lib
|
||||
HSL_LIBS=-lhsl_mc64 -L$(HSL_LIBDIR)
|
||||
HSL_FLAGS= -DHAVE_HSL -I$(HSL_INCDIR)
|
||||
|
||||
# SPRAL package for auction algorithm
|
||||
SPRAL_DIR=/opt/spral/2015.04.20/gnu/6.4.0
|
||||
SPRAL_INCDIR=$(SPRAL_DIR)/include
|
||||
SPRAL_LIBDIR=$(SPRAL_DIR)/lib
|
||||
SPRAL_LIBS=-lspral -L$(SPRAL_LIBDIR)
|
||||
SPRAL_FLAGS=-DHAVE_SPRAL -I$(SPRAL_INCDIR)
|
||||
|
||||
BCM_DIR=/opt/bcm/0.9/gnu/6.4.0
|
||||
BCM_INCDIR=$(BCM_DIR)/include
|
||||
BCM_LIBDIR=$(BCM_DIR)/lib
|
||||
BCM_LDLIBS=-lBCM -L$(BCM_LIBDIR) $(HSL_LIBS) $(SPRAL_LIBS)
|
||||
|
||||
CDEFINES=$(MLDCDEFINES) -I$(BCM_INCDIR)
|
||||
|
||||
LINKOPT=
|
||||
EXEDIR=./runs
|
||||
|
||||
all: mld_d_pde3d
|
||||
|
||||
BCMOBJS= mld_d_bcmatch_aggregator_mod.o mld_d_bcmatch_aggregator_mat_asb.o \
|
||||
mld_d_bcmatch_aggregator_tprol.o mld_d_bcmatch_map_to_tprol.o \
|
||||
mld_daggrmat_unsmth_spmm_asb.o bootCMatch_interface.o
|
||||
|
||||
mld_d_pde3d: mld_d_pde3d.o data_input.o $(BCMOBJS)
|
||||
$(FLINK) $(LINKOPT) mld_d_pde3d.o data_input.o $(BCMOBJS) -o mld_d_pde3d $(MLD_LIBS) $(BCM_LDLIBS) $(PSBLAS_LIBS) $(LDLIBS)
|
||||
/bin/mv mld_d_pde3d $(EXEDIR)
|
||||
|
||||
mld_d_pde3d.o: data_input.o mld_d_bcmatch_aggregator_mod.o
|
||||
mld_d_bcmatch_aggregator_mat_asb.o mld_d_bcmatch_aggregator_tprol.o mld_d_bcmatch_map_to_tprol.o: mld_d_bcmatch_aggregator_mod.o
|
||||
|
||||
check: all
|
||||
cd runs && ./mld_d_pde2d <mld_pde2d.inp && ./mld_s_pde2d<mld_pde2d.inp
|
||||
|
||||
|
||||
clean:
|
||||
/bin/rm -f data_input.o mld_d_pde3d.o *$(.mod) $(BCMOBJS)\
|
||||
$(EXEDIR)/mld_d_pde3d
|
||||
|
||||
verycleanlib:
|
||||
(cd ../..; make veryclean)
|
||||
lib:
|
||||
(cd ../../; make library)
|
||||
|
||||
|
||||
|
||||
@ -1,3 +0,0 @@
|
||||
This directory contains an experimental interface to the BootCMatch
|
||||
aggregator. It may also serve as an example of how to introduce a new
|
||||
aggregation procedure.
|
||||
@ -1,213 +0,0 @@
|
||||
!
|
||||
!
|
||||
! 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
|
||||
use mld_d_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(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_mpk_) :: ctxt, 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_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_
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt,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_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
|
||||
@ -1,400 +0,0 @@
|
||||
!
|
||||
!
|
||||
! MLD2P4 version 2.2
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
!
|
||||
!
|
||||
! The aggregator object hosts the aggregation method for building
|
||||
! the multilevel hierarchy. This variant is based on the hybrid method
|
||||
! presented in
|
||||
!
|
||||
! S. Gratton, P. Henon, P. Jiranek and X. Vasseur:
|
||||
! Reducing complexity of algebraic multigrid by aggregation
|
||||
! Numerical Lin. Algebra with Applications, 2016, 23:501-518
|
||||
!
|
||||
!
|
||||
! sm - class(mld_T_base_smoother_type), allocatable
|
||||
! The current level preconditioner (aka smoother).
|
||||
! parms - type(mld_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 bcm_csr_type_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
|
||||
|
||||
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
|
||||
end module bcm_csr_type_mod
|
||||
|
||||
module mld_d_bcmatch_aggregator_mod
|
||||
use mld_d_base_aggregator_mod
|
||||
use bcm_csr_type_mod
|
||||
|
||||
type, extends(mld_d_base_aggregator_type) :: mld_d_bcmatch_aggregator_type
|
||||
integer(psb_ipk_) :: matching_alg
|
||||
integer(psb_ipk_) :: n_sweeps
|
||||
real(psb_dpk_), allocatable :: w_tmp(:), w_nxt(:)
|
||||
type(bcm_Vector) :: w_par
|
||||
integer(psb_ipk_) :: max_csize
|
||||
integer(psb_ipk_) :: max_nlevels
|
||||
contains
|
||||
procedure, pass(ag) :: bld_tprol => mld_d_bcmatch_aggregator_build_tprol
|
||||
procedure, pass(ag) :: cseti => d_bcmatch_aggr_cseti
|
||||
procedure, pass(ag) :: default => d_bcmatch_aggr_set_default
|
||||
procedure, pass(ag) :: mat_asb => mld_d_bcmatch_aggregator_mat_asb
|
||||
procedure, pass(ag) :: update_next => d_bcmatch_aggregator_update_next
|
||||
procedure, pass(ag) :: bld_wnxt => d_bcmatch_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_bcmatch_aggregator_descr
|
||||
procedure, pass(ag) :: clone => d_bcmatch_aggregator_clone
|
||||
procedure, pass(ag) :: free => d_bcmatch_aggregator_free
|
||||
!!$ procedure, pass(ag) :: default => mld_d_base_aggregator_default
|
||||
procedure, nopass :: fmt => d_bcmatch_aggregator_fmt
|
||||
end type mld_d_bcmatch_aggregator_type
|
||||
|
||||
|
||||
interface
|
||||
subroutine mld_d_bcmatch_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info)
|
||||
import :: mld_d_bcmatch_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, &
|
||||
& psb_ipk_, psb_long_int_k_, mld_dml_parms
|
||||
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(:)
|
||||
type(psb_dspmat_type), intent(out) :: op_prol
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine mld_d_bcmatch_aggregator_build_tprol
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine mld_d_bcmatch_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,ac,&
|
||||
& op_prol,op_restr,info)
|
||||
import :: mld_d_bcmatch_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, &
|
||||
& psb_ipk_, psb_long_int_k_, mld_dml_parms
|
||||
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_), 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
|
||||
end subroutine mld_d_bcmatch_aggregator_mat_asb
|
||||
end interface
|
||||
|
||||
|
||||
interface
|
||||
subroutine mld_d_bcmatch_map_to_tprol(desc_a,ilaggr,nlaggr,valaggr, op_prol,info)
|
||||
import :: mld_d_bcmatch_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, &
|
||||
& psb_ipk_, psb_long_int_k_, mld_dml_parms
|
||||
implicit none
|
||||
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
|
||||
end subroutine mld_d_bcmatch_map_to_tprol
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine mld_daggrmat_unsmth_spmm_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
|
||||
import :: mld_d_bcmatch_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, &
|
||||
& psb_ipk_, psb_long_int_k_, mld_dml_parms
|
||||
implicit none
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
type(psb_desc_type), intent(in) :: desc_a
|
||||
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
|
||||
type(mld_dml_parms), intent(inout) :: parms
|
||||
type(psb_dspmat_type), intent(inout) :: op_prol
|
||||
type(psb_dspmat_type), intent(out) :: ac,op_restr
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine mld_daggrmat_unsmth_spmm_asb
|
||||
end interface
|
||||
|
||||
|
||||
contains
|
||||
|
||||
subroutine d_bld_default_w(ag,nr)
|
||||
use psb_realloc_mod
|
||||
implicit none
|
||||
class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag
|
||||
integer(psb_ipk_), intent(in) :: nr
|
||||
integer(psb_ipk_) :: info
|
||||
call psb_realloc(nr,ag%w_tmp,info)
|
||||
if (info /= psb_success_) return
|
||||
ag%w_tmp = 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(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag
|
||||
|
||||
ag%w_par%size = psb_size(ag%w_tmp)
|
||||
ag%w_par%owns_data = 0
|
||||
if (ag%w_par%size > 0) call set_cloc(ag%w_tmp, ag%w_par)
|
||||
|
||||
end subroutine d_set_default_bcm_w
|
||||
|
||||
subroutine set_cloc(vect,w_par)
|
||||
use iso_c_binding
|
||||
real(psb_dpk_), target :: vect(:)
|
||||
type(bcm_Vector) :: w_par
|
||||
|
||||
w_par%data = c_loc(vect)
|
||||
end subroutine set_cloc
|
||||
|
||||
|
||||
subroutine d_bcmatch_bld_wnxt(ag,ilaggr,valaggr,nx)
|
||||
use psb_realloc_mod
|
||||
implicit none
|
||||
class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag
|
||||
integer(psb_ipk_), intent(in) :: ilaggr(:)
|
||||
real(psb_dpk_), intent(in) :: valaggr(:)
|
||||
integer(psb_ipk_), intent(in) :: nx
|
||||
|
||||
integer(psb_ipk_) :: info,i,j
|
||||
|
||||
call psb_realloc(nx,ag%w_nxt,info)
|
||||
associate(w_nxt => ag%w_nxt, w_tmp=>ag%w_tmp)
|
||||
w_nxt = dzero
|
||||
if (.false.) then
|
||||
do j=1, size(ilaggr)
|
||||
i = ilaggr(j)
|
||||
w_nxt(i) = w_nxt(i) + valaggr(j)*w_tmp(j)
|
||||
end do
|
||||
!write(0,*) 'Old copy ',w_nxt(1:10)
|
||||
else
|
||||
!write(0,*) 'New copy ',nx
|
||||
do i=1, nx
|
||||
w_nxt(i) = w_tmp(i)
|
||||
end do
|
||||
!write(0,*) 'New copy ',w_nxt(1:10)
|
||||
end if
|
||||
end associate
|
||||
|
||||
end subroutine d_bcmatch_bld_wnxt
|
||||
|
||||
function d_bcmatch_aggregator_fmt() result(val)
|
||||
implicit none
|
||||
character(len=32) :: val
|
||||
|
||||
val = "BootCMatch aggregation"
|
||||
end function d_bcmatch_aggregator_fmt
|
||||
|
||||
subroutine d_bcmatch_aggregator_descr(ag,parms,iout,info)
|
||||
implicit none
|
||||
class(mld_d_bcmatch_aggregator_type), intent(in) :: ag
|
||||
type(mld_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_bcmatch_aggregator_descr
|
||||
|
||||
subroutine d_bcmatch_aggregator_update_next(ag,agnext,info)
|
||||
use psb_realloc_mod
|
||||
implicit none
|
||||
class(mld_d_bcmatch_aggregator_type), target, intent(inout) :: ag
|
||||
class(mld_d_base_aggregator_type), target, intent(inout) :: agnext
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
|
||||
!
|
||||
!
|
||||
select type(agnext)
|
||||
class is (mld_d_bcmatch_aggregator_type)
|
||||
agnext%matching_alg = ag%matching_alg
|
||||
agnext%n_sweeps = ag%n_sweeps
|
||||
agnext%max_csize = ag%max_csize
|
||||
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_tmp,info)
|
||||
call agnext%set_c_default_w()
|
||||
class default
|
||||
! What should we do here?
|
||||
end select
|
||||
info = 0
|
||||
end subroutine d_bcmatch_aggregator_update_next
|
||||
|
||||
subroutine d_bcmatch_aggr_cseti(ag,what,val,info,idx)
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_bcmatch_aggregator_type), intent(inout) :: ag
|
||||
character(len=*), intent(in) :: what
|
||||
integer(psb_ipk_), intent(in) :: val
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer, intent(in), optional :: idx
|
||||
integer(psb_ipk_) :: err_act, iwhat
|
||||
character(len=20) :: name='d_bcmatch_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
|
||||
|
||||
end select
|
||||
return
|
||||
end subroutine d_bcmatch_aggr_cseti
|
||||
|
||||
subroutine d_bcmatch_aggr_set_default(ag)
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
class(mld_d_bcmatch_aggregator_type), intent(inout) :: ag
|
||||
character(len=20) :: name='d_bcmatch_aggr_set_default'
|
||||
ag%matching_alg = 0
|
||||
ag%n_sweeps = 1
|
||||
ag%max_nlevels = 36
|
||||
ag%max_csize = 10
|
||||
|
||||
return
|
||||
|
||||
end subroutine d_bcmatch_aggr_set_default
|
||||
|
||||
subroutine d_bcmatch_aggregator_free(ag,info)
|
||||
use iso_c_binding
|
||||
implicit none
|
||||
class(mld_d_bcmatch_aggregator_type), intent(inout) :: ag
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
|
||||
info = 0
|
||||
if (allocated(ag%w_tmp)) deallocate(ag%w_tmp,stat=info)
|
||||
if (info /= 0) return
|
||||
if (allocated(ag%w_nxt)) deallocate(ag%w_nxt,stat=info)
|
||||
if (info /= 0) return
|
||||
ag%w_par%size = 0
|
||||
ag%w_par%data = c_null_ptr
|
||||
ag%w_par%owns_data = 0
|
||||
end subroutine d_bcmatch_aggregator_free
|
||||
|
||||
subroutine d_bcmatch_aggregator_clone(ag,agnext,info)
|
||||
implicit none
|
||||
class(mld_d_bcmatch_aggregator_type), intent(inout) :: ag
|
||||
class(mld_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 (mld_d_bcmatch_aggregator_type)
|
||||
call agnext%set_c_default_w()
|
||||
class default
|
||||
! Should never ever get here
|
||||
info = -1
|
||||
end select
|
||||
end subroutine d_bcmatch_aggregator_clone
|
||||
|
||||
end module mld_d_bcmatch_aggregator_mod
|
||||
@ -1,201 +0,0 @@
|
||||
!
|
||||
!
|
||||
! 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
|
||||
!
|
||||
!
|
||||
|
||||
subroutine mld_d_bcmatch_aggregator_build_tprol(ag,parms,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 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(:)
|
||||
type(psb_dspmat_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_) :: ctxt, 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_csr_type_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_csr_type_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'
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt,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)
|
||||
nr = a%get_nrows()
|
||||
if (psb_size(ag%w_tmp) < 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,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
|
||||
!!$ 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(ctxt,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
|
||||
@ -1,159 +0,0 @@
|
||||
!
|
||||
!
|
||||
! 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
|
||||
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_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_) :: ctxt,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()
|
||||
!
|
||||
ctxt=desc_a%get_context()
|
||||
call psb_info(ctxt,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
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,246 +0,0 @@
|
||||
!
|
||||
!
|
||||
! 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_daggrmat_nosmth_asb.F90
|
||||
!
|
||||
! Subroutine: mld_daggrmat_nosmth_asb
|
||||
! Version: real
|
||||
!
|
||||
! This routine builds a coarse-level matrix A_C from a fine-level matrix A
|
||||
! by using the Galerkin approach, i.e.
|
||||
!
|
||||
! A_C = P_C^T A P_C,
|
||||
!
|
||||
! where P_C is the piecewise constant interpolation operator corresponding
|
||||
! the fine-to-coarse level mapping built by mld_aggrmap_bld.
|
||||
!
|
||||
! The coarse-level matrix A_C is distributed among the parallel processes or
|
||||
! replicated on each of them, according to the value of p%parms%coarse_mat
|
||||
! specified by the user through 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
|
||||
!
|
||||
! For details see
|
||||
! P. D'Ambra, D. di Serafino and S. Filippone, On the development of
|
||||
! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math.,
|
||||
! 57 (2007), 1181-1196.
|
||||
!
|
||||
!
|
||||
! Arguments:
|
||||
! a - type(psb_dspmat_type), input.
|
||||
! The sparse matrix structure containing the local part of
|
||||
! the fine-level matrix.
|
||||
! desc_a - type(psb_desc_type), input.
|
||||
! The communication descriptor of the fine-level matrix.
|
||||
! p - type(mld_d_onelev_type), input/output.
|
||||
! The 'one-level' data structure that will contain the local
|
||||
! part of the matrix to be built as well as the information
|
||||
! concerning the prolongator and its transpose.
|
||||
! parms - type(mld_dml_parms), input
|
||||
! Parameters controlling the choice of algorithm
|
||||
! ac - type(psb_dspmat_type), output
|
||||
! The coarse matrix on output
|
||||
!
|
||||
! ilaggr - integer, dimension(:), input
|
||||
! The mapping between the row indices of the coarse-level
|
||||
! matrix and the row indices of the fine-level matrix.
|
||||
! ilaggr(i)=j means that node i in the adjacency graph
|
||||
! of the fine-level matrix is mapped onto node j in the
|
||||
! adjacency graph of the coarse-level matrix. Note that the indices
|
||||
! are assumed to be shifted so as to make sure the ranges on
|
||||
! the various processes do not overlap.
|
||||
! nlaggr - integer, dimension(:) input
|
||||
! nlaggr(i) contains the aggregates held by process i.
|
||||
! op_prol - type(psb_dspmat_type), input/output
|
||||
! The tentative prolongator on input, the computed prolongator on output
|
||||
!
|
||||
! op_restr - type(psb_dspmat_type), output
|
||||
! The restrictor operator; normally, it is the transpose of the prolongator.
|
||||
!
|
||||
! info - integer, output.
|
||||
! Error code.
|
||||
!
|
||||
!
|
||||
subroutine mld_daggrmat_unsmth_spmm_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
|
||||
use psb_base_mod
|
||||
use mld_d_inner_mod!, mld_protect_name => mld_daggrmat_unsmth_spmm_asb
|
||||
|
||||
implicit none
|
||||
|
||||
! Arguments
|
||||
type(psb_dspmat_type), intent(in) :: a
|
||||
type(psb_desc_type), intent(in) :: desc_a
|
||||
integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:)
|
||||
type(mld_dml_parms), intent(inout) :: parms
|
||||
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
|
||||
integer(psb_ipk_) :: err_act
|
||||
integer(psb_ipk_) :: ctxt,np,me, icomm, ndx, minfo
|
||||
character(len=20) :: name
|
||||
integer(psb_ipk_) :: ierr(5)
|
||||
type(psb_d_coo_sparse_mat) :: ac_coo, tmpcoo
|
||||
type(psb_d_csr_sparse_mat) :: acsr1, acsr2
|
||||
type(psb_dspmat_type) :: am3, am4, tmp_prol
|
||||
integer(psb_ipk_) :: debug_level, debug_unit
|
||||
integer(psb_ipk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
|
||||
& naggr, nzt, naggrm1, naggrp1, i, k
|
||||
|
||||
name='mld_aggrmat_unsmth_spmm_asb'
|
||||
if(psb_get_errstatus().ne.0) return
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
|
||||
ctxt = desc_a%get_context()
|
||||
icomm = desc_a%get_mpic()
|
||||
call psb_info(ctxt, me, np)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
nglob = desc_a%get_global_rows()
|
||||
nrow = desc_a%get_local_rows()
|
||||
ncol = desc_a%get_local_cols()
|
||||
|
||||
|
||||
naggr = nlaggr(me+1)
|
||||
ntaggr = sum(nlaggr)
|
||||
naggrm1 = sum(nlaggr(1:me))
|
||||
naggrp1 = sum(nlaggr(1:me+1))
|
||||
|
||||
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
|
||||
if (info /= psb_success_) goto 9999
|
||||
|
||||
call op_prol%cp_to(acsr1)
|
||||
|
||||
call tmp_prol%mv_from(acsr1)
|
||||
!
|
||||
! Now we have to gather the halo of tmp_prol, and add it to itself
|
||||
! to multiply it by A,
|
||||
!
|
||||
call psb_sphalo(tmp_prol,desc_a,am4,info,&
|
||||
& colcnv=.false.,rowscale=.true.)
|
||||
if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)
|
||||
if (info == psb_success_) call am4%free()
|
||||
if(info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of tmp_prol')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call psb_spspmm(a,tmp_prol,am3,info)
|
||||
if(info /= psb_success_) then
|
||||
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
call tmp_prol%mv_to(tmpcoo)
|
||||
call tmpcoo%transp()
|
||||
|
||||
nzl = tmpcoo%get_nzeros()
|
||||
i=0
|
||||
!
|
||||
! Now we have to fix this. The only rows of B that are correct
|
||||
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
|
||||
!
|
||||
do k=1, nzl
|
||||
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
|
||||
i = i+1
|
||||
tmpcoo%val(i) = tmpcoo%val(k)
|
||||
tmpcoo%ia(i) = tmpcoo%ia(k)
|
||||
tmpcoo%ja(i) = tmpcoo%ja(k)
|
||||
end if
|
||||
end do
|
||||
call tmpcoo%set_nzeros(i)
|
||||
! call tmpcoo%trim()
|
||||
call op_restr%mv_from(tmpcoo)
|
||||
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
|
||||
goto 9999
|
||||
end if
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),&
|
||||
& 'starting sphalo/ rwxtd'
|
||||
|
||||
! op_restr = ((i-wDA)Ptilde)^T
|
||||
call psb_sphalo(am3,desc_a,am4,info,&
|
||||
& colcnv=.false.,rowscale=.true.)
|
||||
if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)
|
||||
if (info == psb_success_) call am4%free()
|
||||
if(info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
! op_restr
|
||||
call psb_sphalo(am3,desc_a,am4,info,&
|
||||
& colcnv=.false.,rowscale=.true.)
|
||||
if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)
|
||||
if (info == psb_success_) call am4%free()
|
||||
if(info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),&
|
||||
& 'starting spspmm 3'
|
||||
call psb_spspmm(op_restr,am3,ac,info)
|
||||
if (info == psb_success_) call am3%free()
|
||||
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),&
|
||||
& 'Done smooth_aggregate '
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
|
||||
return
|
||||
|
||||
end subroutine mld_daggrmat_unsmth_spmm_asb
|
||||
@ -1,72 +0,0 @@
|
||||
|
||||
#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_nlevels,*max_csize , &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_nlevels,max_csize , &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);
|
||||
}
|
||||
@ -1,324 +0,0 @@
|
||||
!
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
module data_input
|
||||
use psb_base_mod, only : psb_spk_, psb_dpk_, psb_ipk_
|
||||
|
||||
interface read_data
|
||||
module procedure read_char, read_int, read_int_array,&
|
||||
& read_double, read_double_array, &
|
||||
& read_single, read_single_array, read_logical,&
|
||||
& string_read_char, string_read_int, string_read_int_array,&
|
||||
& string_read_double, string_read_double_array,&
|
||||
& string_read_single, string_read_single_array, &
|
||||
& string_read_logical
|
||||
end interface read_data
|
||||
interface trim_string
|
||||
module procedure trim_string
|
||||
end interface
|
||||
|
||||
character(len=4096), private :: charbuf
|
||||
character, private, parameter :: def_marker="!"
|
||||
character, private, parameter :: cmt_marker="%"
|
||||
|
||||
contains
|
||||
subroutine get_buffer(file,buffer)
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=*), intent(inout) :: buffer
|
||||
integer :: idx
|
||||
do
|
||||
read(file,'(a)',end=999) buffer
|
||||
buffer = adjustl(buffer)
|
||||
idx=index(charbuf,cmt_marker)
|
||||
if (idx == 1 ) then
|
||||
cycle
|
||||
else
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
999 continue
|
||||
return
|
||||
end subroutine get_buffer
|
||||
|
||||
subroutine read_logical(val,file,marker)
|
||||
logical, intent(out) :: val
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_logical
|
||||
|
||||
subroutine read_char(val,file,marker)
|
||||
character(len=*), intent(out) :: val
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_char
|
||||
|
||||
subroutine read_int(val,file,marker)
|
||||
integer(psb_ipk_), intent(out) :: val
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_int
|
||||
|
||||
subroutine read_int_array(val,file,marker)
|
||||
integer(psb_ipk_), intent(out) :: val(:)
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_int_array
|
||||
|
||||
subroutine read_single(val,file,marker)
|
||||
real(psb_spk_), intent(out) :: val
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_single
|
||||
|
||||
subroutine read_single_array(val,file,marker)
|
||||
real(psb_spk_), intent(out) :: val(:)
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_single_array
|
||||
|
||||
subroutine read_double(val,file,marker)
|
||||
real(psb_dpk_), intent(out) :: val
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_double
|
||||
|
||||
subroutine read_double_array(val,file,marker)
|
||||
real(psb_dpk_), intent(out) :: val(:)
|
||||
integer(psb_ipk_), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
|
||||
call get_buffer(file,charbuf)
|
||||
call read_data(val,charbuf,marker)
|
||||
|
||||
end subroutine read_double_array
|
||||
|
||||
subroutine string_read_char(val,file,marker)
|
||||
character(len=*), intent(out) :: val
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),'(a)') val
|
||||
end subroutine string_read_char
|
||||
|
||||
subroutine string_read_int(val,file,marker)
|
||||
integer(psb_ipk_), intent(out) :: val
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),*) val
|
||||
end subroutine string_read_int
|
||||
|
||||
subroutine string_read_int_array(val,file,marker)
|
||||
integer(psb_ipk_), intent(out) :: val(:)
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),*) val(:)
|
||||
end subroutine string_read_int_array
|
||||
|
||||
subroutine string_read_single(val,file,marker)
|
||||
real(psb_spk_), intent(out) :: val
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),*) val
|
||||
end subroutine string_read_single
|
||||
|
||||
subroutine string_read_single_array(val,file,marker)
|
||||
real(psb_spk_), intent(out) :: val(:)
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),*) val(:)
|
||||
end subroutine string_read_single_array
|
||||
|
||||
subroutine string_read_double(val,file,marker)
|
||||
real(psb_dpk_), intent(out) :: val
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),*) val
|
||||
end subroutine string_read_double
|
||||
|
||||
subroutine string_read_double_array(val,file,marker)
|
||||
real(psb_dpk_), intent(out) :: val(:)
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),*) val(:)
|
||||
end subroutine string_read_double_array
|
||||
|
||||
subroutine string_read_logical(val,file,marker)
|
||||
logical, intent(out) :: val
|
||||
character(len=*), intent(in) :: file
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=1) :: marker_
|
||||
character(len=1024) :: charbuf
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
read(file,'(a)')charbuf
|
||||
charbuf = adjustl(charbuf)
|
||||
idx=index(charbuf,marker_)
|
||||
if (idx == 0) idx = len(charbuf)+1
|
||||
read(charbuf(1:idx-1),*) val
|
||||
end subroutine string_read_logical
|
||||
|
||||
function trim_string(string,marker)
|
||||
character(len=*), intent(in) :: string
|
||||
character(len=1), optional, intent(in) :: marker
|
||||
character(len=len(string)) :: trim_string
|
||||
character(len=1) :: marker_
|
||||
integer(psb_ipk_) :: idx
|
||||
if (present(marker)) then
|
||||
marker_ = marker
|
||||
else
|
||||
marker_ = def_marker
|
||||
end if
|
||||
idx=index(string,marker_)
|
||||
trim_string = adjustl(string(idx:))
|
||||
end function trim_string
|
||||
end module data_input
|
||||
|
||||
@ -1,52 +0,0 @@
|
||||
%%%%%%%%%%% General arguments % Lines starting with % are ignored.
|
||||
CSR ! Storage format CSR COO JAD
|
||||
0100 ! IDIM; domain size. Linear system size is IDIM**2
|
||||
FCG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
|
||||
2 ! ISTOPC
|
||||
00500 ! ITMAX
|
||||
1 ! ITRACE
|
||||
30 ! IRST (restart for RGMRES and BiCGSTABL)
|
||||
1.d-6 ! EPS
|
||||
ML-KCYCLE-FBGS-ILU ! Longer descriptive name for preconditioner (up to 20 chars)
|
||||
ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML
|
||||
%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
|
||||
FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous.
|
||||
1 ! Number of sweeps for smoother
|
||||
0 ! Number of overlap layers for AS preconditioner
|
||||
HALO ! AS restriction operator: NONE HALO
|
||||
NONE ! AS prolongation operator: NONE SUM AVG
|
||||
ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF
|
||||
0 ! Fill level P for ILU(P) and ILU(T,P)
|
||||
1.d-4 ! Threshold T for ILU(T,P)
|
||||
%%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%%
|
||||
NONE ! Second (post) smoother, ignored if NONE
|
||||
1 ! Number of sweeps for (post) smoother
|
||||
0 ! Number of overlap layers for AS preconditioner
|
||||
HALO ! AS restriction operator: NONE HALO
|
||||
NONE ! AS prolongation operator: NONE SUM AVG
|
||||
ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF
|
||||
0 ! Fill level P for ILU(P) and ILU(T,P)
|
||||
1.d-4 ! Threshold T for ILU(T,P)
|
||||
%%%%%%%%%%% Multilevel parameters %%%%%%%%%%%%%%%%
|
||||
KCYCLE ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD
|
||||
1 ! Number of outer sweeps for ML
|
||||
-3 ! Max Number of levels in a multilevel preconditioner; if <0, lib default
|
||||
-3 ! Target coarse matrix size; if <0, lib default
|
||||
UNSMOOTHED ! Type of aggregation: SMOOTHED UNSMOOTHED
|
||||
DEC ! Parallel aggregation: DEC, SYMDEC
|
||||
NATURAL ! Ordering of aggregation NATURAL DEGREE
|
||||
NOFILTER ! Filtering of matrix: FILTER NOFILTER
|
||||
-1.5 ! Coarsening ratio, if < 0 use library default
|
||||
-2 ! Number of thresholds in vector, next line ignored if <= 0
|
||||
0.05 0.025 ! Thresholds
|
||||
-0.0100d0 ! Smoothed aggregation threshold: >= 0.0
|
||||
%%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%%
|
||||
BJAC ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC DEFLT
|
||||
ILU ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS SLU
|
||||
DIST ! Coarsest-level matrix distribution: DIST REPL, DEFLT
|
||||
1 ! Coarsest-level fillin P for ILU(P) and ILU(T,P)
|
||||
1.d-4 ! Coarsest-level threshold T for ILU(T,P)
|
||||
1 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver
|
||||
T ! Use BootCMatch
|
||||
2 ! Matching method: 0 PREIS, 1 MC64, 2 SPRAL (auction)
|
||||
2 ! Pairing sweeps
|
||||
Loading…
Reference in New Issue