From 75d09c6349734d222314296ac497e276a333222c Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 13 Apr 2021 09:29:15 +0200 Subject: [PATCH 1/4] Delete spurious test dirs --- tests/Bcmatch/Makefile | 59 - tests/Bcmatch/README.txt | 3 - .../amg_d_bcmatch_aggregator_mat_asb.f90 | 213 ---- .../Bcmatch/amg_d_bcmatch_aggregator_mod.F90 | 400 ------ .../amg_d_bcmatch_aggregator_tprol.f90 | 201 --- tests/Bcmatch/amg_d_bcmatch_map_to_tprol.f90 | 159 --- tests/Bcmatch/amg_d_pde3d.f90 | 1109 ----------------- .../Bcmatch/amg_daggrmat_unsmth_spmm_asb.f90 | 246 ---- tests/Bcmatch/bootCMatch_interface.c | 72 -- tests/Bcmatch/data_input.f90 | 324 ----- tests/Bcmatch/runs/mld-bcm.inp | 52 - 11 files changed, 2838 deletions(-) delete mode 100644 tests/Bcmatch/Makefile delete mode 100644 tests/Bcmatch/README.txt delete mode 100644 tests/Bcmatch/amg_d_bcmatch_aggregator_mat_asb.f90 delete mode 100644 tests/Bcmatch/amg_d_bcmatch_aggregator_mod.F90 delete mode 100644 tests/Bcmatch/amg_d_bcmatch_aggregator_tprol.f90 delete mode 100644 tests/Bcmatch/amg_d_bcmatch_map_to_tprol.f90 delete mode 100644 tests/Bcmatch/amg_d_pde3d.f90 delete mode 100644 tests/Bcmatch/amg_daggrmat_unsmth_spmm_asb.f90 delete mode 100644 tests/Bcmatch/bootCMatch_interface.c delete mode 100644 tests/Bcmatch/data_input.f90 delete mode 100644 tests/Bcmatch/runs/mld-bcm.inp diff --git a/tests/Bcmatch/Makefile b/tests/Bcmatch/Makefile deleted file mode 100644 index 9c218dd6..00000000 --- a/tests/Bcmatch/Makefile +++ /dev/null @@ -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 =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 diff --git a/tests/Bcmatch/amg_d_bcmatch_aggregator_mod.F90 b/tests/Bcmatch/amg_d_bcmatch_aggregator_mod.F90 deleted file mode 100644 index b5707233..00000000 --- a/tests/Bcmatch/amg_d_bcmatch_aggregator_mod.F90 +++ /dev/null @@ -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 diff --git a/tests/Bcmatch/amg_d_bcmatch_aggregator_tprol.f90 b/tests/Bcmatch/amg_d_bcmatch_aggregator_tprol.f90 deleted file mode 100644 index b3b36d10..00000000 --- a/tests/Bcmatch/amg_d_bcmatch_aggregator_tprol.f90 +++ /dev/null @@ -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 diff --git a/tests/Bcmatch/amg_d_bcmatch_map_to_tprol.f90 b/tests/Bcmatch/amg_d_bcmatch_map_to_tprol.f90 deleted file mode 100644 index 609a29d7..00000000 --- a/tests/Bcmatch/amg_d_bcmatch_map_to_tprol.f90 +++ /dev/null @@ -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 diff --git a/tests/Bcmatch/amg_d_pde3d.f90 b/tests/Bcmatch/amg_d_pde3d.f90 deleted file mode 100644 index ec39a67c..00000000 --- a/tests/Bcmatch/amg_d_pde3d.f90 +++ /dev/null @@ -1,1109 +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. -! -! -! -! File: mld_d_pde3d.f90 -! -! Program: mld_d_pde3d -! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! -! -! The PDE is a general second order equation in 3d -! -! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) -! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f -! dxdx dydy dzdz dx dy dz -! -! with Dirichlet boundary conditions -! u = g -! -! on the unit cube 0<=x,y,z<=1. -! -! -! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. -! -! There are three choices available for data distribution: -! 1. A simple BLOCK distribution -! 2. A ditribution based on arbitrary assignment of indices to processes, -! typically from a graph partitioner -! 3. A 3D distribution in which the unit cube is partitioned -! into subcubes, each one assigned to a process. -! -module mld_d_pde3d_mod - use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& - & psb_dspmat_type, psb_d_vect_type, dzero,& - & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type - - interface - function d_func_3d(x,y,z) result(val) - import :: psb_dpk_ - real(psb_dpk_), intent(in) :: x,y,z - real(psb_dpk_) :: val - end function d_func_3d - end interface - - interface mld_gen_pde3d - module procedure mld_d_gen_pde3d - end interface mld_gen_pde3d - - -contains - - function d_null_func_3d(x,y,z) result(val) - - real(psb_dpk_), intent(in) :: x,y,z - real(psb_dpk_) :: val - - val = dzero - - end function d_null_func_3d - ! - ! functions parametrizing the differential equation - ! - - ! - ! Note: b1, b2 and b3 are the coefficients of the first - ! derivative of the unknown function. The default - ! we apply here is to have them zero, so that the resulting - ! matrix is symmetric/hermitian and suitable for - ! testing with CG and FCG. - ! When testing methods for non-hermitian matrices you can - ! change the B1/B2/B3 functions to e.g. done/sqrt((3*done)) - ! - function b1(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b1 - real(psb_dpk_), intent(in) :: x,y,z - b1=dzero - end function b1 - function b2(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b2 - real(psb_dpk_), intent(in) :: x,y,z - b2=dzero - end function b2 - function b3(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: b3 - real(psb_dpk_), intent(in) :: x,y,z - b3=dzero - end function b3 - function c(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: c - real(psb_dpk_), intent(in) :: x,y,z - c=dzero - end function c - function a1(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a1 - real(psb_dpk_), intent(in) :: x,y,z - a1=done/80 - end function a1 - function a2(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a2 - real(psb_dpk_), intent(in) :: x,y,z - a2=done/80 - end function a2 - function a3(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a3 - real(psb_dpk_), intent(in) :: x,y,z - a3=done/80 - end function a3 - function g(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: g - real(psb_dpk_), intent(in) :: x,y,z - g = dzero - if (x == done) then - g = done - else if (x == dzero) then - g = exp(y**2-z**2) - end if - end function g - - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine mld_d_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,info,& - & f,amold,vmold,imold,partition,nrl,iv) - use psb_base_mod - use psb_util_mod - ! - ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) - ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions - ! u = g - ! - ! on the unit cube 0<=x,y,z<=1. - ! - ! - ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. - ! - implicit none - integer(psb_ipk_) :: idim - type(psb_dspmat_type) :: a - type(psb_d_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - integer(psb_ipk_) :: ctxt, info - character(len=*) :: afmt - procedure(d_func_3d), optional :: f - class(psb_d_base_sparse_mat), optional :: amold - class(psb_d_base_vect_type), optional :: vmold - class(psb_i_base_vect_type), optional :: imold - integer(psb_ipk_), optional :: partition, nrl,iv(:) - - ! Local variables. - - integer(psb_ipk_), parameter :: nb=20 - type(psb_d_csc_sparse_mat) :: acsc - type(psb_d_coo_sparse_mat) :: acoo - type(psb_d_csr_sparse_mat) :: acsr - real(psb_dpk_) :: zt(nb),x,y,z - integer(psb_ipk_) :: m,n,nnz,nr,nt,glob_row,nlr,i,j,ii,ib,k, partition_ - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - ! For 3D partition - integer(psb_ipk_) :: npx,npy,npz, npdims(3),iamx,iamy,iamz,mynx,myny,mynz - integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:) - ! Process grid - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: icoeff - integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:) - real(psb_dpk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_dpk_) :: deltah, sqdeltah, deltah2 - real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb - integer(psb_ipk_) :: err_act - procedure(d_func_3d), pointer :: f_ - character(len=20) :: name, ch_err,tmpfmt - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ctxt, iam, np) - - - if (present(f)) then - f_ => f - else - f_ => d_null_func_3d - end if - - deltah = done/(idim+2) - sqdeltah = deltah*deltah - deltah2 = (2*done)* deltah - - if (present(partition)) then - if ((1<= partition).and.(partition <= 3)) then - partition_ = partition - else - write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' - partition_ = 3 - end if - else - partition_ = 3 - end if - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - - m = idim*idim*idim - n = m - nnz = ((n*7)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - t0 = psb_wtime() - select case(partition_) - case(1) - ! A BLOCK partition - if (present(nrl)) then - nr = nrl - else - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - end if - - nt = nr - call psb_sum(ctxt,nt) - if (nt /= m) then - write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! First example of use of CDALL: specify for each process a number of - ! contiguous rows - ! - call psb_cdall(ctxt,desc_a,info,nl=nr) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(2) - ! A partition defined by the user through IV - - if (present(iv)) then - if (size(iv) /= m) then - write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - else - write(psb_err_unit,*) iam, 'Initialization error: IV not present' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end if - - ! - ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! - call psb_cdall(ctxt,desc_a,info,vg=iv) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(3) - ! A 3-dimensional partition - - ! A nifty MPI function will split the process list - npdims = 0 - call mpi_dims_create(np,3,npdims,info) - npx = npdims(1) - npy = npdims(2) - npz = npdims(3) - - allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) - ! We can reuse idx2ijk for process indices as well. - call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) - ! Now let's split the 3D cube in hexahedra - call dist1Didx(bndx,idim,npx) - mynx = bndx(iamx+1)-bndx(iamx) - call dist1Didx(bndy,idim,npy) - myny = bndy(iamy+1)-bndy(iamy) - call dist1Didx(bndz,idim,npz) - mynz = bndz(iamz+1)-bndz(iamz) - - ! How many indices do I own? - nlr = mynx*myny*mynz - allocate(myidx(nlr)) - ! Now, let's generate the list of indices I own - nr = 0 - do i=bndx(iamx),bndx(iamx+1)-1 - do j=bndy(iamy),bndy(iamy+1)-1 - do k=bndz(iamz),bndz(iamz+1)-1 - nr = nr + 1 - call ijk2idx(myidx(nr),i,j,k,idim,idim,idim) - end do - end do - end do - if (nr /= nlr) then - write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',& - & nr,nlr,mynx,myny,mynz - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - end if - - ! - ! Third example of use of CDALL: specify for each process - ! the set of global indices it owns. - ! - call psb_cdall(ctxt,desc_a,info,vl=myidx) - - case default - write(psb_err_unit,*) iam, 'Initialization error: should not get here' - info = -1 - call psb_barrier(ctxt) - call psb_abort(ctxt) - return - end select - - - if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess - if (info == psb_success_) call psb_geall(xv,desc_a,info) - if (info == psb_success_) call psb_geall(bv,desc_a,info) - - call psb_barrier(ctxt) - talc = psb_wtime()-t0 - - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='allocation rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! - allocate(val(20*nb),irow(20*nb),& - &icol(20*nb),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - - ! loop over rows belonging to current process in a block - ! distribution. - - call psb_barrier(ctxt) - t1 = psb_wtime() - do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) - icoeff = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) - ! x, y, z coordinates - x = (ix-1)*deltah - y = (iy-1)*deltah - z = (iz-1)*deltah - zt(k) = f_(x,y,z) - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 - if (ix == 1) then - zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k) - else - call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y-1,z) - val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 - if (iy == 1) then - zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k) - else - call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y,z-1) - val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 - if (iz == 1) then - zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k) - else - call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - ! term depending on (x,y,z) - val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & - & + c(x,y,z) - call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - ! term depending on (x,y,z+1) - val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 - if (iz == idim) then - zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k) - else - call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y+1,z) - val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 - if (iy == idim) then - zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k) - else - call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x+1,y,z) - val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 - if (ix==idim) then - zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k) - else - call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - end do - call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) - if(info /= psb_success_) exit - zt(:)=dzero - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='insert rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ctxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info,mold=imold) - tcdasb = psb_wtime()-t1 - call psb_barrier(ctxt) - t1 = psb_wtime() - if (info == psb_success_) then - if (present(amold)) then - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) - else - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - end if - end if - call psb_barrier(ctxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) - if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ctxt) - ttot = psb_wtime() - t0 - - call psb_amx(ctxt,talc) - call psb_amx(ctxt,tgen) - call psb_amx(ctxt,tasb) - call psb_amx(ctxt,ttot) - if(iam == psb_root_) then - tmpfmt = a%get_fmt() - write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& - & tmpfmt - write(psb_out_unit,'("-allocation time : ",es12.5)') talc - write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen - write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb - write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb - write(psb_out_unit,'("-total time : ",es12.5)') ttot - - end if - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(ctxt,err_act) - - return - end subroutine mld_d_gen_pde3d - -end module mld_d_pde3d_mod - -program mld_d_pde3d - use psb_base_mod - use mld_prec_mod - use psb_krylov_mod - use psb_util_mod - use data_input - use mld_d_pde3d_mod - use mld_d_bcmatch_aggregator_mod - implicit none - - ! input parameters - character(len=20) :: kmethd, ptype - character(len=5) :: afmt - integer(psb_ipk_) :: idim - - ! miscellaneous - real(psb_dpk_) :: t1, t2, tprec, thier, tslv - - ! sparse matrix and preconditioner - type(psb_dspmat_type) :: a - type(mld_dprec_type) :: prec - type(mld_d_bcmatch_aggregator_type) :: bcmag - ! descriptor - type(psb_desc_type) :: desc_a - ! dense vectors - type(psb_d_vect_type) :: x,b,r - ! parallel environment - integer(psb_ipk_) :: ctxt, iam, np - - ! solver parameters - integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, nlv - integer(psb_long_int_k_) :: amatsize, precsize, descsize - real(psb_dpk_) :: err, resmx, resmxp - - ! Krylov solver data - type solverdata - character(len=40) :: kmethd ! Krylov solver - integer(psb_ipk_) :: istopc ! stopping criterion - integer(psb_ipk_) :: itmax ! maximum number of iterations - integer(psb_ipk_) :: itrace ! tracing - integer(psb_ipk_) :: irst ! restart - real(psb_dpk_) :: eps ! stopping tolerance - end type solverdata - type(solverdata) :: s_choice - - ! preconditioner data - type precdata - - ! preconditioner type - character(len=40) :: descr ! verbose description of the prec - character(len=10) :: ptype ! preconditioner type - - integer(psb_ipk_) :: outer_sweeps ! number of outer sweeps: sweeps for 1-level, - ! AMG cycles for ML - ! general AMG data - character(len=16) :: mlcycle ! AMG cycle type - integer(psb_ipk_) :: maxlevs ! maximum number of levels in AMG preconditioner - - ! AMG aggregation - character(len=16) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED - character(len=16) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC - character(len=16) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE - character(len=16) :: aggr_filter ! filtering: FILTER, NO_FILTER - real(psb_dpk_) :: mncrratio ! minimum aggregation ratio - real(psb_dpk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector - integer(psb_ipk_) :: thrvsz ! size of threshold vector - real(psb_dpk_) :: athres ! smoothed aggregation threshold - integer(psb_ipk_) :: csize ! minimum size of coarsest matrix - logical :: use_bcm ! use BootCMatch - integer(psb_ipk_) :: bcm_alg ! Matching method: 0 PREIS, 1 MC64, 2 SPRAL (auction) - integer(psb_ipk_) :: bcm_sweeps ! Pairing sweeps - ! AMG smoother or pre-smoother; also 1-lev preconditioner - character(len=16) :: smther ! (pre-)smoother type: BJAC, AS - integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps - integer(psb_ipk_) :: novr ! number of overlap layers - character(len=16) :: restr ! restriction over application of AS - character(len=16) :: prol ! prolongation over application of AS - character(len=16) :: solve ! local subsolver type: ILU, MILU, ILUT, - ! UMF, MUMPS, SLU, FWGS, BWGS, JAC - integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization - real(psb_dpk_) :: thr ! threshold for ILUT factorization - - ! AMG post-smoother; ignored by 1-lev preconditioner - character(len=16) :: smther2 ! post-smoother type: BJAC, AS - integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps - integer(psb_ipk_) :: novr2 ! number of overlap layers - character(len=16) :: restr2 ! restriction over application of AS - character(len=16) :: prol2 ! prolongation over application of AS - character(len=16) :: solve2 ! local subsolver type: ILU, MILU, ILUT, - ! UMF, MUMPS, SLU, FWGS, BWGS, JAC - integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization - real(psb_dpk_) :: thr2 ! threshold for ILUT factorization - - ! coarsest-level solver - character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST - character(len=16) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr. - ! mat.); UMF, MUMPS, SLU, ILU, ILUT, MILU - ! (repl. mat.) - character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, - ! MILU, UMF, MUMPS, SLU - integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization - real(psb_dpk_) :: cthres ! threshold for ILUT factorization - integer(psb_ipk_) :: cjswp ! sweeps for GS or JAC coarsest-lev subsolver - - end type precdata - type(precdata) :: p_choice - - ! other variables - integer(psb_ipk_) :: info, i, k - character(len=20) :: name,ch_err - - info=psb_success_ - - - call psb_init(ctxt) - call psb_info(ctxt,iam,np) - - if (iam < 0) then - ! This should not happen, but just in case - call psb_exit(ctxt) - stop - endif - if(psb_get_errstatus() /= 0) goto 9999 - name='mld_d_pde3d' - call psb_set_errverbosity(itwo) - ! - ! Hello world - ! - if (iam == psb_root_) then - write(*,*) 'Welcome to MLD2P4 version: ',mld_version_string_ - write(*,*) 'This is the ',trim(name),' sample program' - end if - - ! - ! get parameters - ! - call get_parms(ctxt,afmt,idim,s_choice,p_choice) - - ! - ! allocate and fill in the coefficient matrix, rhs and initial guess - ! - - call psb_barrier(ctxt) - t1 = psb_wtime() - call mld_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,info) - call psb_barrier(ctxt) - t2 = psb_wtime() - t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='mld_gen_pde3d' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - if (iam == psb_root_) & - & write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 - if (iam == psb_root_) & - & write(psb_out_unit,'(" ")') - ! - ! initialize the preconditioner - ! - call prec%init(ctxt,p_choice%ptype,info) - select case(trim(psb_toupper(p_choice%ptype))) - case ('NONE','NOPREC') - ! Do nothing, keep defaults - - case ('JACOBI','GS','FWGS','FBGS') - ! 1-level sweeps from "outer_sweeps" - call prec%set('smoother_sweeps', p_choice%jsweeps, info) - - case ('BJAC') - call prec%set('smoother_sweeps', p_choice%jsweeps, info) - call prec%set('sub_solve', p_choice%solve, info) - call prec%set('sub_fillin', p_choice%fill, info) - call prec%set('sub_iluthrs', p_choice%thr, info) - - case('AS') - call prec%set('smoother_sweeps', p_choice%jsweeps, info) - call prec%set('sub_ovr', p_choice%novr, info) - call prec%set('sub_restr', p_choice%restr, info) - call prec%set('sub_prol', p_choice%prol, info) - call prec%set('sub_solve', p_choice%solve, info) - call prec%set('sub_fillin', p_choice%fill, info) - call prec%set('sub_iluthrs', p_choice%thr, info) - - case ('ML') - ! multilevel preconditioner - - call prec%set('ml_cycle', p_choice%mlcycle, info) - call prec%set('outer_sweeps', p_choice%outer_sweeps,info) - if (p_choice%csize>0)& - & call prec%set('min_coarse_size', p_choice%csize, info) - if (p_choice%mncrratio>1)& - & call prec%set('min_cr_ratio', p_choice%mncrratio, info) - if (p_choice%maxlevs>0)& - & call prec%set('max_levs', p_choice%maxlevs, info) - if (p_choice%athres >= dzero) & - & call prec%set('aggr_thresh', p_choice%athres, info) - if (p_choice%thrvsz>0) then - do k=1,min(p_choice%thrvsz,size(prec%precv)-1) - call prec%set('aggr_thresh', p_choice%athresv(k), info,ilev=(k+1)) - end do - end if - - call prec%set('aggr_prol', p_choice%aggr_prol, info) - call prec%set('par_aggr_alg', p_choice%par_aggr_alg, info) - call prec%set('aggr_ord', p_choice%aggr_ord, info) - call prec%set('aggr_filter', p_choice%aggr_filter,info) - - - call prec%set('smoother_type', p_choice%smther, info) - call prec%set('smoother_sweeps', p_choice%jsweeps, info) - - select case (psb_toupper(p_choice%smther)) - case ('GS','BWGS','FBGS','JACOBI') - ! do nothing - case default - call prec%set('sub_ovr', p_choice%novr, info) - call prec%set('sub_restr', p_choice%restr, info) - call prec%set('sub_prol', p_choice%prol, info) - call prec%set('sub_solve', p_choice%solve, info) - call prec%set('sub_fillin', p_choice%fill, info) - call prec%set('sub_iluthrs', p_choice%thr, info) - end select - - if (psb_toupper(p_choice%smther2) /= 'NONE') then - call prec%set('smoother_type', p_choice%smther2, info,pos='post') - call prec%set('smoother_sweeps', p_choice%jsweeps2, info,pos='post') - select case (psb_toupper(p_choice%smther2)) - case ('GS','BWGS','FBGS','JACOBI') - ! do nothing - case default - call prec%set('sub_ovr', p_choice%novr2, info,pos='post') - call prec%set('sub_restr', p_choice%restr2, info,pos='post') - call prec%set('sub_prol', p_choice%prol2, info,pos='post') - call prec%set('sub_solve', p_choice%solve2, info,pos='post') - call prec%set('sub_fillin', p_choice%fill2, info,pos='post') - call prec%set('sub_iluthrs', p_choice%thr2, info,pos='post') - end select - end if - - call prec%set('coarse_solve', p_choice%csolve, info) - if (psb_toupper(p_choice%csolve) == 'BJAC') & - & call prec%set('coarse_subsolve', p_choice%csbsolve, info) - call prec%set('coarse_mat', p_choice%cmat, info) - call prec%set('coarse_fillin', p_choice%cfill, info) - call prec%set('coarse_iluthrs', p_choice%cthres, info) - call prec%set('coarse_sweeps', p_choice%cjswp, info) - if (p_choice%use_bcm) then - call prec%set(bcmag,info) - call prec%set('BCM_MATCH_ALG',p_choice%bcm_alg, info) - call prec%set('BCM_SWEEPS',p_choice%bcm_sweeps, info) -!!$ if (p_choice%csize>0) call prec%set('BCM_MAX_CSIZE',p_choice%csize, info) - call prec%set('BCM_MAX_NLEVELS',p_choice%maxlevs, info) - !call prec%set('BCM_W_SIZE',desc_a%get_local_rows(), info,ilev=2) - end if - - end select - - ! build the preconditioner - call psb_barrier(ctxt) - t1 = psb_wtime() - !call psb_set_debug_level(9999) - call prec%hierarchy_build(a,desc_a,info) - !call psb_set_debug_level(0) - thier = psb_wtime()-t1 - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_hierarchy_bld') - goto 9999 - end if - call psb_barrier(ctxt) - t1 = psb_wtime() - call prec%smoothers_build(a,desc_a,info) - tprec = psb_wtime()-t1 - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_smoothers_bld') - goto 9999 - end if - - call psb_amx(ctxt, thier) - call psb_amx(ctxt, tprec) - - if(iam == psb_root_) then - write(psb_out_unit,'(" ")') - write(psb_out_unit,'("Preconditioner: ",a)') trim(p_choice%descr) - write(psb_out_unit,'("Preconditioner time: ",es12.5)')thier+tprec - write(psb_out_unit,'(" ")') - end if - - ! - ! iterative method parameters - ! - call psb_barrier(ctxt) - t1 = psb_wtime() - call psb_krylov(s_choice%kmethd,a,prec,b,x,s_choice%eps,& - & desc_a,info,itmax=s_choice%itmax,iter=iter,err=err,itrace=s_choice%itrace,& - & istop=s_choice%istopc,irst=s_choice%irst) - call psb_barrier(ctxt) - tslv = psb_wtime() - t1 - - call psb_amx(ctxt,tslv) - - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='solver routine' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_barrier(ctxt) - tslv = psb_wtime() - t1 - call psb_amx(ctxt,tslv) - - ! compute residual norms - call psb_geall(r,desc_a,info) - call r%zero() - call psb_geasb(r,desc_a,info) - call psb_geaxpby(done,b,dzero,r,desc_a,info) - call psb_spmm(-done,a,x,done,r,desc_a,info) - resmx = psb_genrm2(r,desc_a,info) - resmxp = psb_geamax(r,desc_a,info) - - amatsize = a%sizeof() - descsize = desc_a%sizeof() - precsize = prec%sizeof() - call psb_sum(ctxt,amatsize) - call psb_sum(ctxt,descsize) - call psb_sum(ctxt,precsize) - call prec%descr(iout=psb_out_unit) - if (iam == psb_root_) then - write(psb_out_unit,'("Computed solution on ",i8," processors")') np - write(psb_out_unit,'("Krylov method : ",a)') trim(s_choice%kmethd) - write(psb_out_unit,'("Preconditioner : ",a)') trim(p_choice%descr) - write(psb_out_unit,'("Iterations to convergence : ",i12)') iter - write(psb_out_unit,'("Relative error estimate on exit : ",es12.5)') err - write(psb_out_unit,'("Number of levels in hierarchy : ",i12)') prec%get_nlevs() - write(psb_out_unit,'("Time to build hierarchy : ",es12.5)') thier - write(psb_out_unit,'("Time to build smoothers : ",es12.5)') tprec - write(psb_out_unit,'("Total time for preconditioner : ",es12.5)') tprec+thier - write(psb_out_unit,'("Time to solve system : ",es12.5)') tslv - write(psb_out_unit,'("Time per iteration : ",es12.5)') tslv/iter - write(psb_out_unit,'("Total time : ",es12.5)') tslv+tprec+thier - write(psb_out_unit,'("Residual 2-norm : ",es12.5)') resmx - write(psb_out_unit,'("Residual inf-norm : ",es12.5)') resmxp - write(psb_out_unit,'("Total memory occupation for A : ",i12)') amatsize - write(psb_out_unit,'("Total memory occupation for DESC_A : ",i12)') descsize - write(psb_out_unit,'("Total memory occupation for PREC : ",i12)') precsize - write(psb_out_unit,'("Storage format for A : ",a )') a%get_fmt() - write(psb_out_unit,'("Storage format for DESC_A : ",a )') desc_a%get_fmt() - - end if - - ! - ! cleanup storage and exit - ! - call psb_gefree(b,desc_a,info) - call psb_gefree(x,desc_a,info) - call psb_spfree(a,desc_a,info) - call prec%free(info) - call psb_cdfree(desc_a,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='free routine' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_exit(ctxt) - stop - -9999 continue - call psb_error(ctxt) - -contains - ! - ! get iteration parameters from standard input - ! - ! - ! get iteration parameters from standard input - ! - subroutine get_parms(icontxt,afmt,idim,solve,prec) - - use psb_base_mod - implicit none - - integer(psb_ipk_) :: icontxt, idim - character(len=*) :: afmt - type(solverdata) :: solve - type(precdata) :: prec - integer(psb_ipk_) :: iam, nm, np, inp_unit - character(len=1024) :: filename - - call psb_info(icontxt,iam,np) - - if (iam == psb_root_) then - if (command_argument_count()>0) then - call get_command_argument(1,filename) - inp_unit = 30 - open(inp_unit,file=filename,action='read',iostat=info) - if (info /= 0) then - write(psb_err_unit,*) 'Could not open file ',filename,' for input' - call psb_abort(icontxt) - stop - else - write(psb_err_unit,*) 'Opened file ',trim(filename),' for input' - end if - else - inp_unit=psb_inp_unit - end if - ! read input data - ! - call read_data(afmt,inp_unit) ! matrix storage format - call read_data(idim,inp_unit) ! Discretization grid size - ! Krylov solver data - call read_data(solve%kmethd,inp_unit) ! Krylov solver - call read_data(solve%istopc,inp_unit) ! stopping criterion - call read_data(solve%itmax,inp_unit) ! max num iterations - call read_data(solve%itrace,inp_unit) ! tracing - call read_data(solve%irst,inp_unit) ! restart - call read_data(solve%eps,inp_unit) ! tolerance - ! preconditioner type - call read_data(prec%descr,inp_unit) ! verbose description of the prec - call read_data(prec%ptype,inp_unit) ! preconditioner type - ! First smoother / 1-lev preconditioner - call read_data(prec%smther,inp_unit) ! smoother type - call read_data(prec%jsweeps,inp_unit) ! (pre-)smoother / 1-lev prec sweeps - call read_data(prec%novr,inp_unit) ! number of overlap layers - call read_data(prec%restr,inp_unit) ! restriction over application of AS - call read_data(prec%prol,inp_unit) ! prolongation over application of AS - call read_data(prec%solve,inp_unit) ! local subsolver - call read_data(prec%fill,inp_unit) ! fill-in for incomplete LU - call read_data(prec%thr,inp_unit) ! threshold for ILUT - ! Second smoother/ AMG post-smoother (if NONE ignored in main) - call read_data(prec%smther2,inp_unit) ! smoother type - call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps - call read_data(prec%novr2,inp_unit) ! number of overlap layers - call read_data(prec%restr2,inp_unit) ! restriction over application of AS - call read_data(prec%prol2,inp_unit) ! prolongation over application of AS - call read_data(prec%solve2,inp_unit) ! local subsolver - call read_data(prec%fill2,inp_unit) ! fill-in for incomplete LU - call read_data(prec%thr2,inp_unit) ! threshold for ILUT - ! general AMG data - call read_data(prec%mlcycle,inp_unit) ! AMG cycle type - call read_data(prec%outer_sweeps,inp_unit) ! number of 1lev/outer sweeps - call read_data(prec%maxlevs,inp_unit) ! max number of levels in AMG prec - call read_data(prec%csize,inp_unit) ! min size coarsest mat - ! aggregation - call read_data(prec%aggr_prol,inp_unit) ! aggregation type - call read_data(prec%par_aggr_alg,inp_unit) ! parallel aggregation alg - call read_data(prec%aggr_ord,inp_unit) ! ordering for aggregation - call read_data(prec%aggr_filter,inp_unit) ! filtering - call read_data(prec%mncrratio,inp_unit) ! minimum aggregation ratio - call read_data(prec%thrvsz,inp_unit) ! size of aggr thresh vector - if (prec%thrvsz > 0) then - call psb_realloc(prec%thrvsz,prec%athresv,info) - call read_data(prec%athresv,inp_unit) ! aggr thresh vector - else - read(inp_unit,*) ! dummy read to skip a record - end if - call read_data(prec%athres,inp_unit) ! smoothed aggr thresh - ! coasest-level solver - call read_data(prec%csolve,inp_unit) ! coarsest-lev solver - call read_data(prec%csbsolve,inp_unit) ! coarsest-lev subsolver - call read_data(prec%cmat,inp_unit) ! coarsest mat layout - call read_data(prec%cfill,inp_unit) ! fill-in for incompl LU - call read_data(prec%cthres,inp_unit) ! Threshold for ILUT - call read_data(prec%cjswp,inp_unit) ! sweeps for GS/JAC subsolver - call read_data(prec%use_bcm,inp_unit) - call read_data(prec%bcm_alg,inp_unit) - call read_data(prec%bcm_sweeps,inp_unit) - if (inp_unit /= psb_inp_unit) then - close(inp_unit) - end if - end if - - call psb_bcast(icontxt,afmt) - call psb_bcast(icontxt,idim) - - call psb_bcast(icontxt,solve%kmethd) - call psb_bcast(icontxt,solve%istopc) - call psb_bcast(icontxt,solve%itmax) - call psb_bcast(icontxt,solve%itrace) - call psb_bcast(icontxt,solve%irst) - call psb_bcast(icontxt,solve%eps) - - call psb_bcast(icontxt,prec%descr) - call psb_bcast(icontxt,prec%ptype) - - ! broadcast first (pre-)smoother / 1-lev prec data - call psb_bcast(icontxt,prec%smther) - call psb_bcast(icontxt,prec%jsweeps) - call psb_bcast(icontxt,prec%novr) - call psb_bcast(icontxt,prec%restr) - call psb_bcast(icontxt,prec%prol) - call psb_bcast(icontxt,prec%solve) - call psb_bcast(icontxt,prec%fill) - call psb_bcast(icontxt,prec%thr) - ! broadcast second (post-)smoother - call psb_bcast(icontxt,prec%smther2) - call psb_bcast(icontxt,prec%jsweeps2) - call psb_bcast(icontxt,prec%novr2) - call psb_bcast(icontxt,prec%restr2) - call psb_bcast(icontxt,prec%prol2) - call psb_bcast(icontxt,prec%solve2) - call psb_bcast(icontxt,prec%fill2) - call psb_bcast(icontxt,prec%thr2) - - ! broadcast AMG parameters - call psb_bcast(icontxt,prec%mlcycle) - call psb_bcast(icontxt,prec%outer_sweeps) - call psb_bcast(icontxt,prec%maxlevs) - - call psb_bcast(icontxt,prec%aggr_prol) - call psb_bcast(icontxt,prec%par_aggr_alg) - call psb_bcast(icontxt,prec%aggr_ord) - call psb_bcast(icontxt,prec%aggr_filter) - call psb_bcast(icontxt,prec%mncrratio) - call psb_bcast(ctxt,prec%thrvsz) - if (prec%thrvsz > 0) then - if (iam /= psb_root_) call psb_realloc(prec%thrvsz,prec%athresv,info) - call psb_bcast(ctxt,prec%athresv) - end if - call psb_bcast(ctxt,prec%athres) - - call psb_bcast(icontxt,prec%csize) - call psb_bcast(icontxt,prec%cmat) - call psb_bcast(icontxt,prec%csolve) - call psb_bcast(icontxt,prec%csbsolve) - call psb_bcast(icontxt,prec%cfill) - call psb_bcast(icontxt,prec%cthres) - call psb_bcast(icontxt,prec%cjswp) - call psb_bcast(ctxt,prec%use_bcm) - call psb_bcast(ctxt,prec%bcm_alg) - call psb_bcast(ctxt,prec%bcm_sweeps) - - end subroutine get_parms - -end program mld_d_pde3d diff --git a/tests/Bcmatch/amg_daggrmat_unsmth_spmm_asb.f90 b/tests/Bcmatch/amg_daggrmat_unsmth_spmm_asb.f90 deleted file mode 100644 index 7bd10a3b..00000000 --- a/tests/Bcmatch/amg_daggrmat_unsmth_spmm_asb.f90 +++ /dev/null @@ -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 diff --git a/tests/Bcmatch/bootCMatch_interface.c b/tests/Bcmatch/bootCMatch_interface.c deleted file mode 100644 index 2b3a5087..00000000 --- a/tests/Bcmatch/bootCMatch_interface.c +++ /dev/null @@ -1,72 +0,0 @@ - -#include -#include - -#include "bcm.h" - -bcm_CSRMatrix bootCMatch(bcm_CSRMatrix *C, int *match_algorithm, int *n_sweeps, int *max_nlevels, int *max_csize, bcm_Vector *w); -bcm_CSRMatrix bootCMatch(bcm_CSRMatrix *C, int *match_algorithm, int *n_sweeps, int *max_nlevels, int *max_csize, bcm_Vector *w){ - bcm_Vector *w_temp; - int info; - //double *w_inp; - //w_inp=bcm_VectorData(w); - - bcm_CSRMatrix *P; - bcm_CSRMatrix *Ac; - int ftcoarse=1; - int cr_it=0, cr_relax_type=0; - double cr_relax_weight=0.0; - // Here I am building Ac but I won't use it. - Ac=bcm_CSRMatchingAgg(C, &w, &P, *match_algorithm, *n_sweeps, *max_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= 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 From 39a9c4e4eda4742d02b61a4e7ce46cb964edb0d2 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 12 May 2021 21:34:59 +0200 Subject: [PATCH 2/4] Fix copyright --- samples/advanced/pdegen/data_input.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/samples/advanced/pdegen/data_input.f90 b/samples/advanced/pdegen/data_input.f90 index b25cdeb0..6b961352 100644 --- a/samples/advanced/pdegen/data_input.f90 +++ b/samples/advanced/pdegen/data_input.f90 @@ -1,14 +1,14 @@ ! ! -! MLD2P4 version 2.2 -! MultiLevel Domain Decomposition Parallel Preconditioners Package -! based on PSBLAS (Parallel Sparse BLAS version 3.5) +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) ! -! (C) Copyright 2008-2018 +! (C) Copyright 2021 ! ! Salvatore Filippone ! Pasqua D'Ambra -! Daniela di Serafino +! Fabio Durastante ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions @@ -18,14 +18,14 @@ ! 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 +! 3. The name of the AMG4PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS ! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR ! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS From 94ca610bff03fe2b32483424fb1ea3c646556e32 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 28 Jun 2021 18:47:35 +0200 Subject: [PATCH 3/4] Do not print matching statistics --- amgprec/amg_d_matchboxp_mod.f90 | 11 +++++++---- amgprec/amg_s_matchboxp_mod.f90 | 11 +++++++---- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/amgprec/amg_d_matchboxp_mod.f90 b/amgprec/amg_d_matchboxp_mod.f90 index ddbcd591..5847e97f 100644 --- a/amgprec/amg_d_matchboxp_mod.f90 +++ b/amgprec/amg_d_matchboxp_mod.f90 @@ -118,6 +118,7 @@ module dmatchboxp_mod module procedure dPMatchBox end interface PMatchBox + logical, parameter, private :: print_statistics=.false. contains subroutine dmatchboxp_build_prol(w,a,desc_a,ilaggr,nlaggr,prol,info,& @@ -421,11 +422,13 @@ contains nlpairs = v(3) end block - if (iam == 0) then - write(0,*) 'Matching statistics: Unmatched nodes ',& - & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs + if (print_statistics) then + if (iam == 0) then + write(0,*) 'Matching statistics: Unmatched nodes ',& + & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs + end if end if - + if (display_out_) then block integer(psb_ipk_) :: idx diff --git a/amgprec/amg_s_matchboxp_mod.f90 b/amgprec/amg_s_matchboxp_mod.f90 index 08c598af..5d3f2266 100644 --- a/amgprec/amg_s_matchboxp_mod.f90 +++ b/amgprec/amg_s_matchboxp_mod.f90 @@ -118,6 +118,7 @@ module smatchboxp_mod module procedure sPMatchBox end interface PMatchBox + logical, parameter, private :: print_statistics=.false. contains subroutine smatchboxp_build_prol(w,a,desc_a,ilaggr,nlaggr,prol,info,& @@ -421,11 +422,13 @@ contains nlpairs = v(3) end block - if (iam == 0) then - write(0,*) 'Matching statistics: Unmatched nodes ',& - & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs + if (print_statistics) then + if (iam == 0) then + write(0,*) 'Matching statistics: Unmatched nodes ',& + & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs + end if end if - + if (display_out_) then block integer(psb_ipk_) :: idx From aba9b297170d09c2f9d98aeb68e8eeb877c65c84 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 15 Jul 2021 11:59:34 +0200 Subject: [PATCH 4/4] Fix samples/simple internal and external docs. --- docs/amg4psblas_1.0-guide.pdf | Bin 1830849 -> 1831453 bytes docs/html/userhtmlsu10.html | 38 +- docs/html/userhtmlsu11.html | 62 +- docs/html/userhtmlsu12.html | 64 +- docs/html/userhtmlsu13.html | 62 +- docs/html/userhtmlsu14.html | 22 +- docs/html/userhtmlsu15.html | 46 +- docs/html/userhtmlsu16.html | 114 ++-- docs/html/userhtmlsu6.html | 32 +- docs/html/userhtmlsu7.html | 34 +- docs/html/userhtmlsu9.html | 720 +++++++++++----------- docs/src/gettingstarted.tex | 32 +- docs/src/userinterface.tex | 17 +- samples/simple/pdegen/amg_dexample_ml.f90 | 19 +- samples/simple/pdegen/amg_sexample_ml.f90 | 19 +- 15 files changed, 650 insertions(+), 631 deletions(-) diff --git a/docs/amg4psblas_1.0-guide.pdf b/docs/amg4psblas_1.0-guide.pdf index 84b4b37aac69e604d7c250b83fecc958ea3ce221..4362844f29fce36f5d3ece4916f81039dd1c8936 100644 GIT binary patch delta 2152 zcmah|Yitx%6z0ypXuH@_=FZOSK)Y3eT}r1jvokw8fGCs}TOJh(q=}fglr7!3uxodz zXkv=VAcZJVG#rhOpuY?%4+EJ1iK0EjkLmk3aj3d> zeXvw*2TNn>FMASl9G)5Cr`v1Jum>u=UBxv;NY(NRRBh%%@JltWv^QO3-}0aqXuQn6 z?`dh84Rs!j;QB~C?Yb*0=ku7lN(%|Fl@Uz)xj)zw$ZyJ$u4(Ytb#}WiA(^HQ zH8d9DxYAVrOxBD&jq^Y4syK<@CLb_0lDba>(x-Xq|!-ZV7FWp;klT-P8_SXJl?pDbC z6!Ph=%u-oayIj937jVlpKz=a9!RfD9f4^x+2}5L3E0Z&3# z{{0w&J=|~*_HE$i!qp=84b=2;74$b?Az$Ps!`kaOxIWg$krg6~CrcnJY{&Y<-K~=i zMOPF(E+^|{b)qa!jFmUe%PpTpL`I^srl_*2$bYBUWe3EsBTZ`DG+NZ8qN`1YChMl$ z++>nh+|;bhlxNLoYMCCH_J1k%gfn8jpLtccJAM*3q9xH0%m}2Mo?Pm_jFoMF_WUgZQ9W$E-L8O9sU$qaS~SX20vD;76fd z6XC5vQD+W50mlc$#~X)ahvq3QQmCO2M)>MthZsxUyD8KnF0dHZ42ffz`f;#jNSwgT ze*iulB4_n{3ulJJ5cBLc_?`OuGB9pf42h-vPL0`g`y2}G6xvYe?rYFa5}1xM=ou!3 zte@e9VX;2=(n)8-G9*ko3jfwKUP@1*&_E%E`0hNsOcJ;(?<6QRY*5e;zu_<(Dv4wJU&%V$ zR1zK_3K(Gqh3Saj+UXGYQ@4piGliLm-+9v^$jveWQO^ji6j~5};6aCYkh!3fT2(0z`Ks`M%JpA%)q`vvSeC)UjiK1Fi}9AKywpg`v| zP=UmCX@{W08mOgEiG)thiH)W%xo>8G4q_mJgudfYdtPi^xS`oe(2%e%;aEBnPSiM7 z0tvsiIhKLKq1PNMiNb|v9Lq%Er$2D4DJXoDzF=lxDiV|1VBdN2Tt2x7$Ipv);G?A@ Ykof&+*m6OPR2aIhmX9Cb+}=|DC$jWQS^xk5 delta 1989 zcmai!Yitx%6vumKyIWYh_QBrWot?&R(ZZ~y)0v&;7DfA@P1W*{heji0OQ&|V-G$vn z+o%y`kWz`3D4fJ75+8{P5kbLGd034`6UC?>_y93Ll&Ez3(4=2f)O%-4jmE?e_m}hd z-*eBocV}PtF8jgx>^2KojODRPEwAOX%B*?Te9LcLWtCe2>uRgQs#p1r4(vHL1at1CeF+(8b_*1-# zw*Yt1#9Q$0tDz9pRYDVoV>ai*_amRA3!;>OC(G~=uPi8<3@!8Uus3N4N)q0EDzg3R zq$a3FGN$MWK~W@o)qHGv_C<00pA9wqeXe~XfJaN=@ZLx*Y^%UY*k6sG|2sX`)(Pke z;q|WWR0!|#8ImB&26Ttj|!97#6X5w1t}u}|?g*sF@XEWqmR!Md=f2njhROOlXCNU^np zvAg;1RAv~)?}&s`o3g1~-pm&o2h0b|fp}v=G5CCT-~lt2r;g5VN@X(Hf}7|u-PqQ= z+3Q9gOc(mySWniq3_cf^4Bp&4klvIkn7#1W_hAfop2vYlQyIETllk6szK~1zY$a(o zFxSOkiX?{I;6Q3PyS0#CN<2Nm7y3<7Sm({YxWs2&a|spdeBLY+(wTeRT|`F><*zB+ zV5+ZR=I}oRGVj)R`-@&6{SBsusf5hjpUT|pl+o$AZqStr@kUt}`8Dam^gt>%FdUb} zBv1R$Ww#H~o!g4>@2zebrxn+nEA3u5>qh^h3;X>GSQ+uCy871}bN7Xc<*??#&^Z|U z0|%gaD0CCd><(SEbH`AKtT0&-vRJaB_KqR$Xx9=|B0m9Hk(CB9u~-xr$Nf#++1_^Y zMlrs9)8~#sX zKmIM(;A5cCe&a`O6ydrOrUa?&6e}pY5UyDZD<-%)My-G?6I>0``4T)$ee*QDIl(no zW!KR?2TLh@D17KMxH7@1%*RK8;eo5 zj^a9s>nW~9p?Cm(E^>>Q+eC;?k`DS#!;O<%n6VDPdg6D#6L2ctOyw;UT@-C7ba|tb z=%ntg6e}riLE+c|C(+H6Fe3!$;2MfmD7@-*_-K+G_4sZmPI3+9d$u{k9>nsS;o2!q z%L<~zNV%SECIEwTJ zVTK60(YE1 zcj#qDK;DyK$@C0MZiHbwP~>D0#->S=Umb=+(_8~{xgX9=b0WiJATUF&)_)C%Gh9v6 z_W88H2%T-DmZAc&>)vx>Rn#TV%8bxyMq-F9JPZ%ba9xi+@uMS95Nn@xEETb5K6ES% zv9nW-rK2dm$*~L+-BRyZNfdo(y<;__=+`BVwG?sM3dg#E9 -

[nextMethod hierarchy_build

-

-

+

call p,info)

-

This method builds the hierarchy of matrices and restriction/prolongation operators for the multilevel preconditioner init and set. -

Arguments

a

type(psb_).

The sparse matrix structure containing the local part of the matrix to be preconditioned. Note that ].

desc_a

type(psb_desc_type).

The communication descriptor of a. See the PSBLAS User’s Guide @@ -136,9 +136,9 @@ class="cmr-12">].

info

integer, ).

Error code. If no error, 0 is returned. See Section 

-

[next] [up]

-

diff --git a/docs/html/userhtmlsu11.html b/docs/html/userhtmlsu11.html index 08501195..da684b86 100644 --- a/docs/html/userhtmlsu11.html +++ b/docs/html/userhtmlsu11.html @@ -10,7 +10,7 @@ -

[nextMethod smoothers_build

-

-

+

call p,imold])

-

This method builds the smoothers and the coarsest-level solvers for the multilevel preconditioner  5.3). -

Arguments

@@ -262,7 +262,7 @@ class="td11">

a

type(psb_).

The sparse matrix structure containing the local part of the matrix to be preconditioned. Note that ].

desc_a

type(psb_desc_type).

The communication descriptor of a. See the PSBLAS User’s Guide @@ -150,9 +150,9 @@ class="cmr-12">].

info

integer, ).

Error code. If no error, 0 is returned. See Section 7 for details.

amold

class(psb_optional.

The desired dynamic type for internal matrix components; this allows e.g. running on GPUs; it needs not be the same on all @@ -205,9 +205,9 @@ class="cmr-12">].

vmold

class(psb_optional.

The desired dynamic type for internal vector components; this allows e.g. running on GPUs.

imold

class(psb_i_base_vect_typeoptional.

The desired dynamic type for internal integer vector components; this allows e.g. running on GPUs.

-

[next] [up]

-

diff --git a/docs/html/userhtmlsu12.html b/docs/html/userhtmlsu12.html index 32e905c0..4f4c7dec 100644 --- a/docs/html/userhtmlsu12.html +++ b/docs/html/userhtmlsu12.html @@ -10,7 +10,7 @@ -

[next5.5 Method build

-

-

+

call p,imold])

-

This method builds the preconditioner p according to the requirements made by the @@ -83,7 +83,7 @@ class="cmtt-12">smoothers_build, whose nomenclature would however be somewhat unnatural when dealing with simple one-level preconditioners. -

Arguments

a

type(psb_).

The sparse matrix structure containing the local part of the matrix to be preconditioned. Note that ].

desc_a

type(psb_desc_type).

The communication descriptor of a. See the PSBLAS User’s Guide @@ -158,9 +158,9 @@ class="cmr-12">].

info

integer, ).

Error code. If no error, 0 is returned. See Section 7 for details.

amold

class(psb_optional.

The desired dynamic type for internal matrix components; this allows e.g. running on GPUs; it needs not be the same on all @@ -213,9 +213,9 @@ class="cmr-12">].

vmold

class(psb_optional.

The desired dynamic type for internal vector components; this allows e.g. running on GPUs.

imold

class(psb_i_base_vect_typeoptional.

The desired dynamic type for internal integer vector components; this allows e.g. running on GPUs.

-

The method can be used to build multilevel preconditioners too. @@ -272,7 +272,7 @@ class="cmr-12">The method can be used to build multilevel preconditioners too.

[next] [up]

-

diff --git a/docs/html/userhtmlsu13.html b/docs/html/userhtmlsu13.html index 6f194dd5..aec509c3 100644 --- a/docs/html/userhtmlsu13.html +++ b/docs/html/userhtmlsu13.html @@ -10,7 +10,7 @@ -

[next5.6 Method apply

-

-

+

call p,work])

-

This method computes y = psb_krylov and hence it is completely transparent to the user. -

Arguments

x

type(kind_parameter), dimension(:), intent(in)—.

The local part of the vector x. Note that , single/double class="cmr-12">precision version of AMG4PSBLAS under use.

y

type(kind_parameter), dimension(:), intent(out)—.

The local part of the vector y. Note that , single/double class="cmr-12">precision version of AMG4PSBLAS under use.

desc_a

type(psb_desc_type).

The communication descriptor associated to the matrix to be preconditioned.

info

integer, ).

Error code. If no error, 0 is returned. See Section 7 for details.

trans

character(lenin).

If trans = 1).

work

type(kind_parameter), dimension(:), optional, target—.

Workspace. Its size should be at least 4 

-

[next] [up]

-

diff --git a/docs/html/userhtmlsu14.html b/docs/html/userhtmlsu14.html index 9dab129c..d2f9804d 100644 --- a/docs/html/userhtmlsu14.html +++ b/docs/html/userhtmlsu14.html @@ -10,7 +10,7 @@ -

[next5.7 Method free

-

-

+

call p,info)

-

This method deallocates the preconditioner data structure p. -

Arguments

info

integer, ).

Error code. If no error, 0 is returned. See Section 

-

[next] [up]

-

diff --git a/docs/html/userhtmlsu15.html b/docs/html/userhtmlsu15.html index 8c3cfa76..07ff43a3 100644 --- a/docs/html/userhtmlsu15.html +++ b/docs/html/userhtmlsu15.html @@ -10,7 +10,7 @@ -

[next5.8 Method descr

-

-

+

call p verbosity])

-

This method prints a description of the preconditioner p to the standard output or to a @@ -68,7 +68,7 @@ class="cmtt-12">build, have been called. -

Arguments

info

integer, ).

Error code. If no error, 0 is returned. See Section 7 for details.

iout

integer, optional.

The id of the file where the preconditioner description will be printed; the default is the standard output.

root

integer, optional.

The id of the process where the preconditioner description will be printed; the default is psb_root_.

verbosity

integer, optional.

The verbosity level of the description. Default value is 0. For values higher than 0, it prints out further information, e.g., for @@ -174,7 +174,7 @@ class="td11">

-

[next] [up]

-

diff --git a/docs/html/userhtmlsu16.html b/docs/html/userhtmlsu16.html index 6b6ab93f..122093bf 100644 --- a/docs/html/userhtmlsu16.html +++ b/docs/html/userhtmlsu16.html @@ -10,7 +10,7 @@ -

[prev]

class="cmr-12">5.9 Auxiliary Methods -

Various functionalities are implemented as additional methods of the preconditioner object. -

+

5.9.1 Method: dump
-

-

+

call p,global_num])

-

Dump on file. -

Arguments

info

integer, ).

Error code. If no error, 0 is returned. See Section 7 for details.

amold

class(psb_optional.

The desired dynamic type for internal matrix components; this allows e.g. running on GPUs; it needs not be the same on all @@ -134,15 +134,15 @@ class="cmr-12">.

-

+

5.9.2 Method: clone
-

-

+

call p)
-

Create a (deep) copy of the preconditioner object. -

Arguments

pout

type(amg_).

The copy of the preconditioner data structure. Note that x must @@ -199,9 +199,9 @@ class="cmr-12">, single/double precision class="cmr-12">version of AMG4PSBLAS under use.

info

integer, ).

Error code. If no error, 0 is returned. See Section for details.

-

+

5.9.3 Method: sizeof
-

-

+

sz =

global

logical, optional.

Whether the global or local preconditioner memory occupation is desired. Default: .

Return memory footprint in bytes. -

+

5.9.4 Method: allocate_wrk
-

-

+

call p vmold])

-

Allocate internal work vectors. Each application of the preconditioner uses a number of work vectors which are allocated internally as necessary; therefore allocation and @@ -326,7 +326,7 @@ class="cmr-12">argument to the apply method. -

Arguments

info

integer, ).

Error code. If no error, 0 is returned. See Section 7 for details.

vmold

class(psb_optional.

The desired dynamic type for internal vector components; this allows e.g. running on GPUs.

-

+

5.9.5 Method: free_wrk
-

-

+

call p(info)

-

Deallocate internal work vectors. -

Arguments

info

integer, ).

Error code. If no error, 0 is returned. See Section . ! with basic smoothed aggregation, 1 hybrid forward/backward ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level ! solver -  call P%init(’ML’,info) +  call P%init(ctxt,’ML’,info) ! ! build the preconditioner   call P%hierarchy_build(A,desc_A,info) @@ -169,8 +169,8 @@ class="cmr-12">. ! set the solver parameters and the initial guess   ... ... ! -! solve Ax=b with preconditioned CG -  call psb_krylov(’CG’,A,P,b,x,tol,desc_A,info) +! solve Ax=b with preconditioned FCG +  call psb_krylov(’FCG’,A,P,b,x,tol,desc_A,info)   ... ... ! ! deallocate the preconditioner @@ -296,7 +296,7 @@ class="cmr-12">. ! build a V-cycle preconditioner with 1 block-Jacobi sweep (with ! ILU(0) on the blocks) as pre- and post-smoother, and 8  block-Jacobi ! sweeps (with ILU(0) on the blocks) as coarsest-level solver -  call P%init(’ML’,info) +  call P%init(ctxt,’ML’,info)   call P%set(’SMOOTHER_TYPE’,’BJAC’,info)   call P%set(’COARSE_SOLVE’,’BJAC’,info)   call P%set(’COARSE_SWEEPS’,8,info) @@ -323,25 +323,27 @@ class="content">setup of a multilevel preconditioner based on the default decoup

-

+

 ... ...
 ! build a W-cycle preconditioner with 2 hybrid Gauss-Seidel sweeps
 ! as pre- and post-smoother, a distributed coarsest
 ! matrix, and MUMPS as coarsest-level solver
-  call P%init(’ML’,info)
+  call P%init(ctxt,’ML’,info)
   call P%set(’PAR_AGGR_ALG’,’COUPLED’,info)
-call P%set(’AGGR_TYPE’,’MATCHBOXP’,info)
-call P%set(’AGGR_SIZE’,8,info)
+  call P%set(’AGGR_TYPE’,’MATCHBOXP’,info)
+  call P%set(’AGGR_SIZE’,8,info)
   call P%set(’ML_CYCLE’,’WCYCLE’,info)
   call P%set(’SMOOTHER_TYPE’,’FBGS’,info)
   call P%set(’SMOOTHER_SWEEPS’,2,info)
   call P%set(’COARSE_SOLVE’,’KRM’,info)
+  call P%set(’COARSE_MAT’,’DIST’,info)
+  call P%set(’KRM_METHOD’,’FCG’,info)
   call P%hierarchy_build(A,desc_A,info)
   call P%smoothers_build(A,desc_A,info)
 ... ...
-

+


Listing 3: setup of a multilevel preconditioner based on the coupled coarsening using @@ -353,7 +355,7 @@ weighted matching
-

@@ -361,18 +363,18 @@ weighted matching
-

+

 ... ...
 ! set RAS with overlap 2 and ILU(0) on the local blocks
-  call P%init(’AS’,info)
+  call P%init(ctxt,’AS’,info)
   call P%set(’SUB_OVR’,2,info)
   call P%bld(A,desc_A,info)
 ... ...
 ! solve Ax=b with preconditioned BiCGSTAB
   call psb_krylov(’BICGSTAB’,A,P,b,x,tol,desc_A,info)
-

+


Listing 4: setup of a one-level Schwarz preconditioner.
@@ -383,7 +385,7 @@ class="content">setup of a one-level Schwarz preconditioner.

[next] [up]

-

diff --git a/docs/html/userhtmlsu7.html b/docs/html/userhtmlsu7.html index f6633c56..57eae31c 100644 --- a/docs/html/userhtmlsu7.html +++ b/docs/html/userhtmlsu7.html @@ -10,7 +10,7 @@ -

[prev]

class="cmr-12">4.2 GPU example -

The code discussed here shows how to set up a program exploiting the combined GPU capabilities of PSBLAS and AMG4PSBLAS. The code example is available in the @@ -36,14 +36,14 @@ class="cmr-12">capabilities of PSBLAS and AMG4PSBLAS. The code example is availa class="cmr-12">source distribution directory amg4psblas/examples/gpu. -

First of all, we need to include the appropriate modules and declare some auxiliary variables: -

@@ -51,7 +51,7 @@ class="cmr-12">variables:
-

+

 program amg_dexample_gpu
   use psb_base_mod
@@ -69,7 +69,7 @@ program amg_dexample_gpu
   type(psb_i_vect_gpu)       :: igmold
  
-

+


Listing 5: setup of a GPU-enabled test program part one.
@@ -77,7 +77,7 @@ class="content">setup of a GPU-enabled test program part one.

In this particular example we are choosing to employ a HLG data structure for @@ -85,7 +85,7 @@ class="cmr-12">data structure for class="cmr-12">sparse matrices on GPUs; for more information please refer to the PSBLAS-EXT users’ guide. -

We then have to initialize the GPU environment, and pass the appropriate MOLD variables to the build methods (see also the PSBLAS and PSBLAS-EXT users’ @@ -94,7 +94,7 @@ class="cmr-12">guides). -

@@ -102,7 +102,7 @@ class="cmr-12">guides).
-

+

   call psb_init(ctxt)
   call psb_info(ctxt,iam,np)
@@ -116,7 +116,7 @@ class="cmr-12">guides).
   call prec%smoothers_build(a,desc_a,info, amold=agmold, vmold=vgmold, imold=igmold)
  
-

+


Listing 6: setup of a GPU-enabled test program part two.
@@ -124,7 +124,7 @@ class="content">setup of a GPU-enabled test program part two.

Finally, we convert the input matrix, the descriptor and the vectors to use a GPU-enabled internal storage format. We then preallocate the preconditioner @@ -135,7 +135,7 @@ class="cmr-12">GPU environment -

@@ -143,7 +143,7 @@ class="cmr-12">GPU environment
-

+

   call desc_a%cnv(mold=igmold)
   call a%cscnv(info,mold=agmold)
@@ -167,7 +167,7 @@ class="cmr-12">GPU environment
   stop
  
-

+


Listing 7: setup of a GPU-enabled test program part three.
@@ -175,7 +175,7 @@ class="content">setup of a GPU-enabled test program part three.

It is very important to employ solvers that are suited to the GPU, i.e. solvers that do NOT employ triangular system solve kernels. Solvers that satisfy this constraint @@ -193,7 +193,7 @@ class="cmtt-12">INVT

  • AINV
    • -

    and their ℓ1 ). class="newline" />

    Remark 3. Many of the coarsest-level solvers cannot be used with both replicated +class="cmr-12">Many of the coarsest-level solvers apply to a specific coarsest-matrix and distributed coarsest-matrix layouts; therefore, setting the solver after the layout +class="cmr-12">layout; therefore, setting the solver after the layout may change the layout to either may change the layout. Similarly, setting the layout after the solver may change the +class="cmr-12">distributed or replicated. Similarly, setting the layout after the solver may change the solver. -

    More precisely, UMFPACK and SuperLU require the coarsest-level matrix to be replicated, while SuperLUsingle Gauss-Seidel sweep, respectively, which are very poor solvers. -

    On the other hand, the distributed layout can be used with any solver but UMFPACK and SuperLU; therefore, if any of these two solvers has already been @@ -514,13 +514,13 @@ class="cmr-12">_Dist or KRM have been class="cmr-12">previously set, the coarsest-level solver is changed to the default sequential solver. -

    In a parallel setting with many cores, we suggest to the users to change the default coarsest solver for using the KRM choice, i.e. a parallel distributed iterative solution of the coarsest system based on Krylov methods. -

    Remark 4. The argument idx . -

    @@ -553,7 +553,7 @@ class="cmr-12">.
    -

    +

    +class="td11">

    1

    Parameters defining the aggregation algorithm.

    @@ -954,7 +954,7 @@ class="content">Parameters defining the aggregation algorithm.

    +

    what

    		type

    val

    deault

    cots

    ’ML_CYCLE’

    character(len=*)

    ’VCYCLE’ -

    ’WCYCLE’ -

    ’KCYCLE’ -

    ’ADD’

    ’VCYCLE’

    Multilevel cycle: V-cycle, W-cycle, K-cycle, +class="td11">

    Multilevel cycle: V-cycle, W-cycle, K-cycle, and additive composition.

    ’CYCLE_SWEEPS’

    integer

    Any integer -

    number

    Any integer +

    number ≥ 1

    1

    Number of multilevel cycles.

    Number of multilevel cycles.

    -

    @@ -669,7 +669,7 @@ applied.
    -

    +

    @@ -683,7 +683,7 @@ id="TBL-5-4">

    what

    		type

    val

    deault

    cots

    ’MIN_COARSE_SIZE_PER_PROCESS’

    integer

    Any number -

    Any number +

    > 0

    200

    Coarse size threshold per process. The +class="td11">

    200

    Coarse size threshold per process. The aggregation stops if the global number of variables of the computed coarsest matrix is lower than or equal to this threshold @@ -737,17 +737,17 @@ multiplied by the number of processes

    ’MIN_COARSE_SIZE’

    integer

    Any number -

    Any number +

    > 0

    -1

    Coarse size threshold. The aggregation +class="td11">

    -1

    Coarse size threshold. The aggregation stops if the global number of variables of the computed coarsest matrix is lower than or equal @@ -760,17 +760,17 @@ class="cmtt-10x-x-109">’.

    ’MIN_CR_RATIO’

    real

    Any number -

    Any number +

    > 1

    1.5

    Minimum coarsening +class="td11">

    1.5

    Minimum coarsening ratio. The aggregation stops if the ratio between the global matrix dimensions at two consecutive levels is lower than or @@ -778,41 +778,41 @@ equal to this threshold (see Note).

    ’MAX_LEVS’

    integer

    Any integer -

    number

    Any integer +

    number > 1

    20

    Maximum number of levels. The +class="td11">

    20

    Maximum number of levels. The aggregation stops if the number of levels reaches this value (see Note).

    ’PAR_AGGR_ALG’

    character(len=*)

    ’DEC’, ’SYMDEC’, ’COUPLED’

    ’DEC’

    Parallel aggregation algorithm. -

    the

    Parallel aggregation algorithm. +

    the SYMDEC option applies decoupled aggregation to the sparsity pattern of T .

    ’AGGR_TYPE’

    character(len=*)

    ’SOC1’, @@ -843,11 +843,11 @@ class="cmtt-10x-x-109">’, class="cmtt-10x-x-109">’MATCHBOXP’

    ’SOC1’

    Type of aggregation algorithm: currently, +class="td11">

    Type of aggregation algorithm: currently, for the decoupled aggregation we implement two measures of strength of connection, the @@ -863,22 +863,22 @@ href="userhtmlli5.html#XMatchBoxP">9].

    ’AGGR_SIZE’

    integer

    Any integer -

    power of 2, +class="td11">

    Any integer +

    power of 2, with aggr_size ≥ 2

    4

    Maximum size of aggregates when the +class="td11">

    4

    Maximum size of aggregates when the coupled aggregation based on matching is applied. For aggressive coarsening with size of aggregate larger than 8 @@ -890,16 +890,16 @@ class="cmtt-10x-x-109">’MATCHBOXP’

    ’AGGR_PROL’

    character(len=*)

    ’SMOOTHED’, @@ -907,11 +907,11 @@ class="cmtt-10x-x-109">’, class="cmtt-10x-x-109">’UNSMOOTHED’

    ’SMOOTHED’

    Prolongator used by the aggregation +class="td11">

    Prolongator used by the aggregation algorithm: smoothed or unsmoothed (i.e., tentative prolongator).

    Parameters defining the aggregation algorithm (continued).

    @@ -1097,7 +1097,7 @@ class="content">Parameters defining the aggregation algorithm (continued). -

    +

    @@ -1111,7 +1111,7 @@ id="TBL-7-4">

    @@ -1478,7 +1482,7 @@ class="content">Parameters defining the solver at the coarsest level (continued) -

    @@ -1486,7 +1490,7 @@ class="content">Parameters defining the solver at the coarsest level (continued)
    -

    +

    @@ -1500,7 +1504,7 @@ id="TBL-8-4">

    +class="td11">

    The stopping tolerance.

    Additional parameters defining the solver at the coarsest level. -

    @@ -1940,7 +1944,7 @@ class="content">Additional parameters defining the solver at the coarsest level.
    -

    +

    what

    		type

    val

    deault

    cots

    ’AGGR_ORD’

    character(len=*)

    ’NATURAL’ -

    ’DEGREE’

    ’NATURAL’

    Initial ordering of indices for +class="td11">

    Initial ordering of indices for the decoupled aggregation algorithm: either natural ordering or sorted by descending degrees of the nodes in the @@ -1023,44 +1023,44 @@ matrix graph.

    ’AGGR_THRESH’

    real(kind_parameter)

    Any real -

    number 

    Any real +

    number ∈ [0,1]

    0.01

    The threshold

    0.01

    The threshold θ in the strength of connection algorithm. See also the note at the bottom of this table.

    ’AGGR_FILTER’

    character(len=*)

    ’FILTER’ -

    ’NOFILTER’

    ’NOFILTER’

    Matrix used in computing the smoothed +class="td11">

    Matrix used in computing the smoothed prolongator: filtered or unfiltered.

    what

    		type

    val

    deault

    cots

    ’COARSE_MAT’

    character(len=*)

    ’DIST’ -

    ’REPL’

    ’REPL’

    Coarsest matrix layout: distributed among the +class="td11">

    Coarsest matrix layout: distributed among the processes or replicated on each of them.

    ’COARSE_SOLVE’

    character(len=*)

    ’MUMPS’ -

    ’UMF’ -

    ’SLU’ -

    ’SLUDIST’ -

    ’ILU’ +

    ’JACOBI’ -

    ’GS’ -

    ’BJAC’ -

    ’KRM’ -

    ’L1-JACOBI’ -

    ’L1-BJAC’ -

    ’L1-FBGS’

    See Note.

    Solver used at the coarsest level: sequential +class="td11">

    See Note.

    Solver used at the coarsest level: sequential LU from MUMPS, UMFPACK, or SuperLU (plus triangular solve); distributed LU from MUMPS or SuperLU_Dist (plus triangular solve); @@ -1262,53 +1266,53 @@ double precision.

    ’COARSE_SUBSOLVE’

    character(len=*)

    ’ILU’ -

    ’ILUT’ -

    ’MILU’ -

    ’MUMPS’ -

    ’SLU’ -

    ’UMF’ -

    ’INVT’ -

    ’INVK’ -

    ’AINV’

    See Note.

    Solver for the diagonal blocks of the coarsest +class="td11">

    See Note.

    Solver for the diagonal blocks of the coarsest matrix, in case the block Jacobi solver is chosen as coarsest-level solver: ILU(p), ILU(ILU otherwise.

    what

    		type

    val

    deault

    cos

    ’COARSE_SWEEPS’

    integer

    Any +class="td11">

    Any integer -

    number

    number > 0

    10

    Number of sweeps when

    10

    Number of sweeps when JACOBI, GS or BJAC is @@ -1413,39 +1417,39 @@ chosen as coarsest-level solver.

    ’COARSE_FILLIN’

    integer

    Any +class="td11">

    Any integer -

    number

    number ≥ 0

    0

    Fill-in level

    0

    Fill-in level p of the ILU factorizations and first fill-in for the approximate inverses.

    ’COARSE_ILUTHRS’

    real(kind_parameter)

    Any real -

    number

    Any real +

    number ≥ 0

    0

    Drop tolerance

    0

    Drop tolerance t in the ILU(p,t) factorization and first drop-tolerance for the approximate inverses.

    what

    		type

    val

    deault

    cots

    ’BJAC_STOP’

    character(len=*)

    ’FALSE’ -

    ’TRUE’

    ’FALSE’

    Select whether to use a stopping criterion for +class="td11">

    Select whether to use a stopping criterion for the Block-Jacobi method used as a coarse solver.

    ’BJAC_TRACE’

    character(len=*)

    ’FALSE’ -

    ’TRUE’

    ’FALSE’

    Select whether to print a trace for the +class="td11">

    Select whether to print a trace for the calculated residual for the Block-Jacobi method used as a coarse solver.

    ’BJAC_ITRACE’

    integer

    Any integer -

    Any integer +

    > 0

    -1

    Number of iterations after which a trace is to +class="td11">

    -1

    Number of iterations after which a trace is to be printed.

    ’BJAC_RESCHECK’

    integer

    Any integer -

    Any integer +

    > 0

    -1

    Number of iterations after which a residual is +class="td11">

    -1

    Number of iterations after which a residual is to be calculated.

    ’BJAC_STOPTOL’

    real(kind_parameter)

    Any real -

    Any real +

    < 1

    0

    Tolerance for the stopping criterion on the +class="td11">

    0

    Tolerance for the stopping criterion on the residual.

    ’KRM_METHOD’

    character(len=*)

    ’CG’ -

    ’FCG’ -

    ’CGS’ -

    ’CGR’ -

    ’BICG’ -

    ’BICGSTAB’ -

    ’BICGSTABL’ -

    ’RGMRES’

    ’FCG’

    A string that defines the iterative method to +class="td11">

    A string that defines the iterative method to be used when employing a Krylov method ’

    ’KRM_KPREC’

    character(len=*)

    Table 

    Table 1

    ’BJAC’

    The one-level +class="td11">

    The one-level preconditioners from the Table 1 can be used for the coarse Krylov solver.

    ’KRM_SUB_SOLVE’

    character(len=*)

    Table 

    Table 5

    ’ILU’

    Solver for the diagonal blocks of the coarsest +class="td11">

    Solver for the diagonal blocks of the coarsest matrix preconditioner, in case the block Jacobi solver is chosen as 5 applies here.

    ’KRM_GLOBAL’

    character(len=*)

    ’TRUE’, @@ -1795,46 +1799,46 @@ class="cmtt-10x-x-109">’, class="cmtt-10x-x-109">’FALSE’

    ’FALSE’

    Choose between a global Krylov solver, all +class="td11">

    Choose between a global Krylov solver, all unknowns on a single node, or a distributed one. The default choice is the distributed solver.

    ’KRM_EPS’

    real(kind_parameter)

    Real

    Real < 1

    10

    10-6

    The stopping tolerance.

    ’KRM_IRST’

    integer

    Integer -

    Integer +

    ≥ 1

    30

    An integer specifying the restart parameter. +class="td11">

    30

    An integer specifying the restart parameter. This is employed for the BiCGSTABL or RGMRES @@ -1842,16 +1846,16 @@ methods, otherwise it is ignored.

    ’KRM_ISTOPC’

    integer

    Integers +class="td11">

    Integers 1,2,3

    2

    If

    2

    If 1 then the method uses the normwise backward error in the infinity norm; if 2, the @@ -1864,32 +1868,32 @@ for the details.

    ’KRM_ITMAX’

    integer

    Integer -

    Integer +

    ≥ 1

    40

    The maximum number of iterations to +class="td11">

    40

    The maximum number of iterations to perform.

    ’KRM_ITRACE’

    integer

    Integer -

    Integer +

    ≥ 0

    -1

    If

    -1

    If > 0 print out an informational message about convergence every

    ’KRM_FILLIN’

    integer

    Integer -

    Integer +

    ≥ 0

    0

    Fill-in level

    0

    Fill-in level p of the ILU factorizations and first fill-in for the approximate inverses.

    @@ -2274,7 +2278,7 @@ preconditioner. -

    @@ -2282,7 +2286,7 @@ preconditioner.
    -

    +

    what

    		type

    val

    deault

    cots

    ’SMOOTHER_TYPE’

    character(len=*)

    ’JACOBI’ -

    ’GS’ -

    ’BGS’ -

    ’BJAC’ -

    ’AS’ -

    ’L1-JACOBI’ -

    ’L1-BJAC’ -

    ’L1-FBGS’

    ’FBGS’

    Type of smoother used in the multilevel preconditioner: point-Jacobi, hybrid @@ -2056,21 +2060,21 @@ class="cmr-10">–hybrid (forward) Gauss-Seidel, class="cmmi-10">ℓ1-point-Jacobi and Additive Schwarz. -

    It is ignored by one-level preconditioners.

    ’SUB_SOLVE’

    character(len=*)

    ’JACOBI’ @@ -2078,47 +2082,47 @@ class="cmtt-10">’ class="cmtt-10">’GS’ -

    ’BGS’ -

    ’ILU’ -

    ’ILUT’ -

    ’MILU’ -

    ’MUMPS’ -

    ’SLU’ -

    ’UMF’ -

    ’INVT’ -

    ’INVK’ -

    ’AINV’

    GS and BGS multilevel class="cmr-10">preconditioners, respectively -

    ILU for block-Jacobi and Additive Schwarz class="cmr-10">one-level preconditioners

    The local solver to be used with the smoother or one-level preconditioner (see @@ -2196,22 +2200,22 @@ class="cmr-10">Note for details on hybrid Gauss-Seidel.

    ’SMOOTHER_SWEEPS’

    integer

    Any integer -

    number ≥ 0

    1

    Number of sweeps of the smoother or one-level preconditioner. In the multilevel @@ -2238,22 +2242,22 @@ class="cmr-10">respectively.

    ’SUB_OVR’

    integer

    Any integer -

    number ≥ 0

    1

    Number of overlap layers, for Additive Schwarz only.

    what

    		type

    val

    deault

    cots

    ’SUB_RESTR’

    character(len=*)

    ’HALO’ -

    ’NONE’

    ’HALO’

    Type of restriction operator, for Additive Schwarz only: ’NONE’ for neglecting it. -

    Note that HALO must be chosen for the @@ -2372,28 +2376,28 @@ class="cmr-10">its RAS variant.

    ’SUB_PROL’

    character(len=*)

    ’SUM’ -

    ’NONE’

    ’NONE’

    Type of prolongation operator, for Additive Schwarz only: ’ for neglecting them. -

    Note that ’SUMfor its RAS variant.

    ’SUB_FILLIN’

    integer

    Any integer -

    number ≥ 0

    0

    Fill-in level p of the incomplete LU @@ -2449,25 +2453,25 @@ class="cmr-10">factorizations.

    ’SUB_ILUTHRS’

    real(kind_parameter)

    Any real number ≥ 0

    0

    Drop tolerance t in the ILU(factorization.

    ’MUMPS_LOC_GLOB’

    character(len=*)

    ’LOCAL_SOLVER’ -

    ’GLOBAL_SOLVER’

    ’GLOBAL_SOLVER’

    Whether MUMPS should be used as a distributed solver, or as a serial solver acting @@ -2510,19 +2514,19 @@ class="cmr-10">process.

    ’MUMPS_IPAR_ENTRY’

    integer

    Any integer number

    0

    Set an entry in the MUMPS integer control array, as chosen via the argument.

    ’MUMPS_RPAR_ENTRY’

    real

    Any real number

    0

    Set an entry in the MUMPS real control array, as chosen via the Parameters defining the smoother or the details of the one-level -

    [next] [up]

    -

    diff --git a/docs/src/gettingstarted.tex b/docs/src/gettingstarted.tex index 5dba3495..04b07490 100644 --- a/docs/src/gettingstarted.tex +++ b/docs/src/gettingstarted.tex @@ -174,7 +174,7 @@ the corresponding codes are available in \verb|samples/simple/fileread/|. ! with basic smoothed aggregation, 1 hybrid forward/backward ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level ! solver - call P%init('ML',info) + call P%init(ctxt,'ML',info) ! ! build the preconditioner call P%hierarchy_build(A,desc_A,info) @@ -184,8 +184,8 @@ the corresponding codes are available in \verb|samples/simple/fileread/|. ! set the solver parameters and the initial guess ... ... ! -! solve Ax=b with preconditioned CG - call psb_krylov('CG',A,P,b,x,tol,desc_A,info) +! solve Ax=b with preconditioned FCG + call psb_krylov('FCG',A,P,b,x,tol,desc_A,info) ... ... ! ! deallocate the preconditioner @@ -229,7 +229,7 @@ stop ! with basic smoothed aggregation, 1 hybrid forward/backward ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level ! solver - call P%init('ML',info) + call P%init(ctxt,'ML',info) ! ! build the preconditioner call P%hierarchy_build(A,desc_A,info) @@ -239,8 +239,8 @@ stop ! set the solver parameters and the initial guess ... ... ! -! solve Ax=b with preconditioned CG - call psb_krylov('CG',A,P,b,x,tol,desc_A,info) +! solve Ax=b with preconditioned FCG + call psb_krylov('FCG',A,P,b,x,tol,desc_A,info) ... ... ! ! deallocate the preconditioner @@ -307,7 +307,7 @@ boundary conditions are also available in the directory \verb|samples/simple/pde ! build a V-cycle preconditioner with 1 block-Jacobi sweep (with ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi ! sweeps (with ILU(0) on the blocks) as coarsest-level solver -call P%init('ML',info) +call P%init(ctxt,'ML',info) call P%set('SMOOTHER_TYPE','BJAC',info) call P%set('COARSE_SOLVE','BJAC',info) call P%set('COARSE_SWEEPS',8,info) @@ -323,7 +323,7 @@ call P%smoothers_build(A,desc_A,info) ! build a V-cycle preconditioner with 1 block-Jacobi sweep (with ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi ! sweeps (with ILU(0) on the blocks) as coarsest-level solver - call P%init('ML',info) + call P%init(ctxt,'ML',info) call P%set('SMOOTHER_TYPE','BJAC',info) call P%set('COARSE_SOLVE','BJAC',info) call P%set('COARSE_SWEEPS',8,info) @@ -345,7 +345,7 @@ call P%smoothers_build(A,desc_A,info) !2 hybrid Gauss-Seidel sweeps as pre- and post-smoother, !and parallel flexible Conjugate Gradient coupled with the block-Jacobi !preconditioner having ILU(0) on the blocks as coarsest solver. -call P%init('ML',info) +call P%init(ctxt,'ML',info) call P%set('PAR_AGGR_ALG','COUPLED',info) call P%set('AGGR_TYPE','MATCHBOXP',info) call P%set('AGGR_SIZE',8,info) @@ -353,6 +353,8 @@ call P%set('ML_CYCLE','WCYCLE',info) call P%set('SMOOTHER_TYPE','FBGS',info) call P%set('SMOOTHER_SWEEPS',2,info) call P%set('COARSE_SOLVE','KRM',info) +call P%set('COARSE_MAT','DIST',info) +call P%set('KRM_METHOD','FCG',info) call P%hierarchy_build(A,desc_A,info) call P%smoothers_build(A,desc_A,info) \end{minted} @@ -365,14 +367,16 @@ call P%smoothers_build(A,desc_A,info) ! build a W-cycle preconditioner with 2 hybrid Gauss-Seidel sweeps ! as pre- and post-smoother, a distributed coarsest ! matrix, and MUMPS as coarsest-level solver - call P%init('ML',info) + call P%init(ctxt,'ML',info) call P%set('PAR_AGGR_ALG','COUPLED',info) -call P%set('AGGR_TYPE','MATCHBOXP',info) -call P%set('AGGR_SIZE',8,info) + call P%set('AGGR_TYPE','MATCHBOXP',info) + call P%set('AGGR_SIZE',8,info) call P%set('ML_CYCLE','WCYCLE',info) call P%set('SMOOTHER_TYPE','FBGS',info) call P%set('SMOOTHER_SWEEPS',2,info) call P%set('COARSE_SOLVE','KRM',info) + call P%set('COARSE_MAT','DIST',info) + call P%set('KRM_METHOD','FCG',info) call P%hierarchy_build(A,desc_A,info) call P%smoothers_build(A,desc_A,info) ... ... @@ -387,7 +391,7 @@ call P%set('AGGR_SIZE',8,info) \ifpdf \begin{minted}[breaklines=true,bgcolor=bg,fontsize=\small]{fortran} ! build a one-level RAS with overlap 2 and ILU(0) on the local blocks. -call P%init('AS',info) +call P%init(ctxt,'AS',info) call P%set('SUB_OVR',2,info) call P%build(A,desc_A,info) ... ... @@ -401,7 +405,7 @@ call P%build(A,desc_A,info) \begin{verbatim} ... ... ! set RAS with overlap 2 and ILU(0) on the local blocks - call P%init('AS',info) + call P%init(ctxt,'AS',info) call P%set('SUB_OVR',2,info) call P%bld(A,desc_A,info) ... ... diff --git a/docs/src/userinterface.tex b/docs/src/userinterface.tex index 2a496079..007a3469 100644 --- a/docs/src/userinterface.tex +++ b/docs/src/userinterface.tex @@ -176,9 +176,10 @@ Similar considerations apply to the point-Jacobi, Gauss-Seidel and block-Jacobi coarsest-level solvers, and shortcuts are available in this case too (see Table~\ref{tab:p_coarse_1}). \\ -\textbf{Remark 3.} Many of the coarsest-level solvers cannot be used -with both replicated and distributed coarsest-matrix layouts; -therefore, setting the solver after the layout may change the layout. +\textbf{Remark 3.} Many of the coarsest-level solvers apply to a +specific coarsest-matrix layout; +therefore, setting the solver after the layout may change the layout +to either distributed or replicated. Similarly, setting the layout after the solver may change the solver. More precisely, UMFPACK and SuperLU require the coarsest-level @@ -399,8 +400,14 @@ the parameter \texttt{ilev}.} \\ \fortinline|'COARSE_SOLVE'| & \fortinline|character(len=*)| & \fortinline|'MUMPS'| \par \fortinline|'UMF'| \par \fortinline|'SLU'| \par \fortinline|'SLUDIST'| \par - \fortinline|'JACOBI'| \par \fortinline|'GS'| \par \fortinline|'BJAC'| \par \fortinline|'KRM'| - \fortinline|'L1-JACOBI'| \par \fortinline|'L1-BJAC'| \par \fortinline|'L1-FBGS'| + \fortinline|'ILU'| \par + \fortinline|'JACOBI'| \par + \fortinline|'GS'| \par + \fortinline|'BJAC'| \par + \fortinline|'KRM'| \par + \fortinline|'L1-JACOBI'| \par + \fortinline|'L1-BJAC'| \par + \fortinline|'L1-FBGS'| & See~Note. & Solver used at the coarsest level: sequential LU from MUMPS, UMFPACK, or SuperLU diff --git a/samples/simple/pdegen/amg_dexample_ml.f90 b/samples/simple/pdegen/amg_dexample_ml.f90 index 64bdc4fc..c9c18fbd 100644 --- a/samples/simple/pdegen/amg_dexample_ml.f90 +++ b/samples/simple/pdegen/amg_dexample_ml.f90 @@ -38,7 +38,7 @@ ! File: amg_dexample_ml.f90 ! ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. The solver is CG, coupled with one of the +! PDE with Dirichlet BCs. The solver is FCG, coupled with one of the ! following multi-level preconditioner, as explained in Section 4.1 of ! the AMG4PSBLAS User's and Reference Guide: ! @@ -121,7 +121,8 @@ program amg_dexample_ml real(psb_dpk_) :: resmx, resmxp real(psb_dpk_) :: t1, t2, tprec character(len=5) :: afmt='CSR' - character(len=20) :: name, kmethod + character(len=20) :: name + character(len=20), parameter :: kmethod='FCG' ! initialize the parallel environment @@ -177,7 +178,6 @@ program amg_dexample_ml ! solver call P%init(ctxt,'ML',info) - kmethod = 'CG' case(2) @@ -188,15 +188,16 @@ program amg_dexample_ml call P%init(ctxt,'ML',info) call P%set('SMOOTHER_TYPE','BJAC',info) call P%set('COARSE_SOLVE','BJAC',info) + call P%set('COARSE_SUBSOLVE','ILU',info) call P%set('COARSE_SWEEPS',8,info) - kmethod = 'CG' case(3) - ! initialize a W-cycle preconditioner based on the coupled aggregation relying on matching, - ! with maximum size of aggregates equal to 8 and smoothed prolongators, - ! 2 hybrid forward/backward GS sweeps as pre/post-smoother, a distributed coarsest - ! matrix, and preconditioned Flexible Conjugate Gradient as coarsest-level solver + ! initialize a W-cycle preconditioner based on the coupled aggregation + ! relying on matching, with maximum size of aggregates equal to 8 + ! and smoothed prolongators, 2 hybrid forward/backward GS sweeps + ! as pre/post-smoother, a distributed coarsest matrix, + ! and preconditioned Flexible Conjugate Gradient as coarsest-level solver call P%init(ctxt,'ML',info) call P%set('PAR_AGGR_ALG','COUPLED',info) @@ -206,7 +207,7 @@ program amg_dexample_ml call P%set('SMOOTHER_SWEEPS',2,info) call P%set('COARSE_SOLVE','KRM',info) call P%set('COARSE_MAT','DIST',info) - kmethod = 'CG' + call P%set('KRM_METHOD','FCG',info) end select diff --git a/samples/simple/pdegen/amg_sexample_ml.f90 b/samples/simple/pdegen/amg_sexample_ml.f90 index 0f2734ea..e3d023e1 100644 --- a/samples/simple/pdegen/amg_sexample_ml.f90 +++ b/samples/simple/pdegen/amg_sexample_ml.f90 @@ -38,7 +38,7 @@ ! File: amg_sexample_ml.f90 ! ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. The solver is CG, coupled with one of the +! PDE with Dirichlet BCs. The solver is FCG, coupled with one of the ! following multi-level preconditioner, as explained in Section 4.1 of ! the AMG4PSBLAS User's and Reference Guide: ! @@ -121,7 +121,8 @@ program amg_sexample_ml real(psb_spk_) :: resmx, resmxp real(psb_dpk_) :: t1, t2, tprec character(len=5) :: afmt='CSR' - character(len=20) :: name, kmethod + character(len=20) :: name + character(len=20), parameter :: kmethod='FCG' ! initialize the parallel environment @@ -177,7 +178,6 @@ program amg_sexample_ml ! solver call P%init(ctxt,'ML',info) - kmethod = 'CG' case(2) @@ -188,15 +188,16 @@ program amg_sexample_ml call P%init(ctxt,'ML',info) call P%set('SMOOTHER_TYPE','BJAC',info) call P%set('COARSE_SOLVE','BJAC',info) + call P%set('COARSE_SUBSOLVE','ILU',info) call P%set('COARSE_SWEEPS',8,info) - kmethod = 'CG' case(3) - ! initialize a W-cycle preconditioner based on the coupled aggregation relying on matching, - ! with maximum size of aggregates equal to 8 and smoothed prolongators, - ! 2 hybrid forward/backward GS sweeps as pre/post-smoother, a distributed coarsest - ! matrix, and preconditioned Flexible Conjugate Gradient as coarsest-level solver + ! initialize a W-cycle preconditioner based on the coupled aggregation + ! relying on matching, with maximum size of aggregates equal to 8 + ! and smoothed prolongators, 2 hybrid forward/backward GS sweeps + ! as pre/post-smoother, a distributed coarsest matrix, + ! and preconditioned Flexible Conjugate Gradient as coarsest-level solver call P%init(ctxt,'ML',info) call P%set('PAR_AGGR_ALG','COUPLED',info) @@ -206,7 +207,7 @@ program amg_sexample_ml call P%set('SMOOTHER_SWEEPS',2,info) call P%set('COARSE_SOLVE','KRM',info) call P%set('COARSE_MAT','DIST',info) - kmethod = 'CG' + call P%set('KRM_METHOD','FCG',info) end select