From ea27db555037f767084eb10dc80e4fc311cce44c Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 16 Mar 2018 20:06:44 +0100 Subject: [PATCH] make Utilities::MPI::sum() support SparseMatrix --- include/deal.II/base/mpi.h | 16 +++- include/deal.II/base/mpi.templates.h | 28 +++++- include/deal.II/lac/sparse_matrix.h | 46 ++++++++- source/base/mpi.inst.in | 6 ++ tests/mpi/collective_sparse_matrix_01.cc | 95 +++++++++++++++++++ ...ollective_sparse_matrix_01.mpirun=1.output | 2 + ...ollective_sparse_matrix_01.mpirun=3.output | 2 + 7 files changed, 188 insertions(+), 7 deletions(-) create mode 100644 tests/mpi/collective_sparse_matrix_01.cc create mode 100644 tests/mpi/collective_sparse_matrix_01.mpirun=1.output create mode 100644 tests/mpi/collective_sparse_matrix_01.mpirun=3.output diff --git a/include/deal.II/base/mpi.h b/include/deal.II/base/mpi.h index 2dd20725c3..7084595734 100644 --- a/include/deal.II/base/mpi.h +++ b/include/deal.II/base/mpi.h @@ -52,7 +52,7 @@ DEAL_II_NAMESPACE_OPEN //Forward type declarations to allow MPI sums over tensorial types template class Tensor; template class SymmetricTensor; - +template class SparseMatrix; namespace Utilities { @@ -196,6 +196,20 @@ namespace Utilities sum (const Tensor &local, const MPI_Comm &mpi_communicator); + /** + * Perform an MPI sum of the entries of a SparseMatrix. + * + * @note @p local and @p global should have the same sparsity + * pattern and it should be the same for all MPI processes. + * + * @relatesalso SparseMatrix + */ + template + void + sum (const SparseMatrix &local, + const MPI_Comm &mpi_communicator, + SparseMatrix &global); + /** * Return the maximum over all processors of the value @p t. This function * is collective over all processors given in the diff --git a/include/deal.II/base/mpi.templates.h b/include/deal.II/base/mpi.templates.h index 5c8cc75a96..8302d7cdc9 100644 --- a/include/deal.II/base/mpi.templates.h +++ b/include/deal.II/base/mpi.templates.h @@ -24,6 +24,7 @@ #include #include #include +#include #include @@ -235,7 +236,7 @@ namespace Utilities const MPI_Comm &mpi_communicator) { Tensor sums; - dealii::Utilities::MPI::sum(local, mpi_communicator, sums); + sum(local, mpi_communicator, sums); return sums; } @@ -253,7 +254,7 @@ namespace Utilities entries[i] = local[ local.unrolled_to_component_indices(i) ]; Number global_entries[ SymmetricTensor::n_independent_components ]; - dealii::Utilities::MPI::sum( entries, mpi_communicator, global_entries ); + sum( entries, mpi_communicator, global_entries ); SymmetricTensor global; for (unsigned int i=0; i< n_entries; ++i) @@ -264,6 +265,29 @@ namespace Utilities + template + void + sum (const SparseMatrix &local, + const MPI_Comm &mpi_communicator, + SparseMatrix &global) + { + Assert(local.get_sparsity_pattern() == global.get_sparsity_pattern(), + ExcMessage("The sparsity pattern of the local and the global matrices should match.")); +#ifdef DEAL_II_WITH_MPI + // makes use of the fact that the matrix stores its data in a + // contiguous array. + sum(ArrayView (&local.val[0], local.n_nonzero_elements()), + mpi_communicator, + ArrayView(&global.val[0], global.n_nonzero_elements())); +#else + (void)mpi_communicator; + if (!PointerComparison::equal (&local, &global)) + global = local; +#endif + } + + + template T max (const T &t, const MPI_Comm &mpi_communicator) diff --git a/include/deal.II/lac/sparse_matrix.h b/include/deal.II/lac/sparse_matrix.h index 7c87cd7022..42edaa941f 100644 --- a/include/deal.II/lac/sparse_matrix.h +++ b/include/deal.II/lac/sparse_matrix.h @@ -24,6 +24,9 @@ #include #include #include +#ifdef DEAL_II_WITH_MPI +#include +#endif #include @@ -34,6 +37,15 @@ template class Vector; template class FullMatrix; template class BlockMatrixBase; template class SparseILU; +#ifdef DEAL_II_WITH_MPI +namespace Utilities +{ + namespace MPI + { + template void sum (const SparseMatrix &, const MPI_Comm &, SparseMatrix &); + } +} +#endif #ifdef DEAL_II_WITH_TRILINOS namespace TrilinosWrappers @@ -974,8 +986,14 @@ public: * classes instead, since they are tailored better to a sparse matrix * structure. */ - number operator () (const size_type i, - const size_type j) const; + const number &operator () (const size_type i, + const size_type j) const; + + /** + * In contrast to the one above, this function allows modifying the object. + */ + number &operator () (const size_type i, + const size_type j); /** * This function is mostly like operator()() in that it returns the value of @@ -1628,6 +1646,13 @@ private: */ template friend class SparseMatrixIterators::Iterator; template friend class SparseMatrixIterators::Accessor; + +#ifdef DEAL_II_WITH_MPI + /** + * Give access to internal datastructures to perform MPI operations. + */ + template friend void Utilities::MPI::sum (const SparseMatrix &, const MPI_Comm &, SparseMatrix &); +#endif }; #ifndef DOXYGEN @@ -1868,8 +1893,21 @@ SparseMatrix::operator /= (const number factor) template inline -number SparseMatrix::operator () (const size_type i, - const size_type j) const +const number &SparseMatrix::operator () (const size_type i, + const size_type j) const +{ + Assert (cols != nullptr, ExcNotInitialized()); + Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry, + ExcInvalidIndex(i,j)); + return val[cols->operator()(i,j)]; +} + + + +template +inline +number &SparseMatrix::operator () (const size_type i, + const size_type j) { Assert (cols != nullptr, ExcNotInitialized()); Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry, diff --git a/source/base/mpi.inst.in b/source/base/mpi.inst.in index a20a7fbcaf..704c9c780d 100644 --- a/source/base/mpi.inst.in +++ b/source/base/mpi.inst.in @@ -14,6 +14,12 @@ // --------------------------------------------------------------------- +for (S : REAL_SCALARS) +{ + template + void sum (const SparseMatrix &, const MPI_Comm &, SparseMatrix &); +} + for (S : MPI_SCALARS) { template diff --git a/tests/mpi/collective_sparse_matrix_01.cc b/tests/mpi/collective_sparse_matrix_01.cc new file mode 100644 index 0000000000..6845ff91e6 --- /dev/null +++ b/tests/mpi/collective_sparse_matrix_01.cc @@ -0,0 +1,95 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test Utilities::MPI::sum() for sparse matrix + +#include "../tests.h" +#include "../testmatrix.h" +#include +#include + + + +void test() +{ + const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + const unsigned int size = 50; + + FDMatrix testproblem (size, size); + unsigned int dim = (size-1) * (size-1); + + SparsityPattern sparsity(dim, dim, size); + testproblem.five_point_structure(sparsity); + sparsity.compress(); + + SparseMatrix matrix(sparsity); + const double val = std::pow(10,myid); + for (SparsityPattern::const_iterator it = sparsity.begin(); + it != sparsity.end(); ++it) + { + const auto i = (*it).row(); + const auto j = (*it).column(); + matrix.add(i, j, -val); + matrix.add(i, i, val); + } + + // compare with FullMatrix: + FullMatrix full(dim,dim); + full.copy_from(matrix); + + // deallog << "Local:" << std::endl; + // matrix.print(deallog.get_file_stream()); + Utilities::MPI::sum(matrix, MPI_COMM_WORLD, matrix); + // deallog << "Global:" << std::endl; + // matrix.print(deallog.get_file_stream()); + + Utilities::MPI::sum(full, MPI_COMM_WORLD, full); + + for (SparsityPattern::const_iterator it = sparsity.begin(); + it != sparsity.end(); ++it) + { + const auto i = (*it).row(); + const auto j = (*it).column(); + AssertThrow (matrix(i,j) == full(i,j), + ExcMessage(std::to_string(matrix(i,j)) + + " != " + + std::to_string(full(i,j)) + + " for i=" + std::to_string(i) + + " j=" + std::to_string(j) + )); + } + + deallog << "Ok" << std::endl; +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); + + // MPILogInitAll log; + // test(); + + if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) + { + initlog(); + test(); + } + else + test(); + + return 0; +} diff --git a/tests/mpi/collective_sparse_matrix_01.mpirun=1.output b/tests/mpi/collective_sparse_matrix_01.mpirun=1.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/mpi/collective_sparse_matrix_01.mpirun=1.output @@ -0,0 +1,2 @@ + +DEAL::Ok diff --git a/tests/mpi/collective_sparse_matrix_01.mpirun=3.output b/tests/mpi/collective_sparse_matrix_01.mpirun=3.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/mpi/collective_sparse_matrix_01.mpirun=3.output @@ -0,0 +1,2 @@ + +DEAL::Ok -- 2.39.5