From 6be05fa06147de6b65241c105a2d10418b15a80c Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 19 Mar 2019 13:39:57 +0100 Subject: [PATCH] add SparsityTools::gather_sparsity_pattern() --- doc/news/changes/minor/20190319DenisDavydov_b | 3 + .../deal.II/lac/dynamic_sparsity_pattern.h | 6 + include/deal.II/lac/sparsity_tools.h | 28 +++++ source/lac/dynamic_sparsity_pattern.cc | 19 +++ source/lac/sparsity_tools.cc | 111 ++++++++++++++++++ tests/sparsity/dynamic_sparsity_pattern_19.cc | 105 +++++++++++++++++ ...ynamic_sparsity_pattern_19.mpirun=3.output | 71 +++++++++++ 7 files changed, 343 insertions(+) create mode 100644 doc/news/changes/minor/20190319DenisDavydov_b create mode 100644 tests/sparsity/dynamic_sparsity_pattern_19.cc create mode 100644 tests/sparsity/dynamic_sparsity_pattern_19.mpirun=3.output diff --git a/doc/news/changes/minor/20190319DenisDavydov_b b/doc/news/changes/minor/20190319DenisDavydov_b new file mode 100644 index 0000000000..7f728e91ef --- /dev/null +++ b/doc/news/changes/minor/20190319DenisDavydov_b @@ -0,0 +1,3 @@ +New: Add SparsityTools::gather_sparsity_pattern(). +
+(Denis Davydov, 2019/03/19) diff --git a/include/deal.II/lac/dynamic_sparsity_pattern.h b/include/deal.II/lac/dynamic_sparsity_pattern.h index ed2a200179..372f7fadf5 100644 --- a/include/deal.II/lac/dynamic_sparsity_pattern.h +++ b/include/deal.II/lac/dynamic_sparsity_pattern.h @@ -518,6 +518,12 @@ public: size_type row_length(const size_type row) const; + /** + * Clear all entries stored in a specific row. + */ + void + clear_row(const size_type row); + /** * Access to column number field. Return the column number of the @p * indexth entry in @p row. diff --git a/include/deal.II/lac/sparsity_tools.h b/include/deal.II/lac/sparsity_tools.h index e428bef3ff..d681b84796 100644 --- a/include/deal.II/lac/sparsity_tools.h +++ b/include/deal.II/lac/sparsity_tools.h @@ -292,6 +292,34 @@ namespace SparsityTools const MPI_Comm & mpi_comm, const IndexSet & myrange); + /** + * Gather rows in a dynamic sparsity pattern over MPI. + * The function is similar to SparsityTools::distribute(), + * however instead of distributing sparsity stored in non-owned + * rows on this MPI process, this function will gather sparsity + * from other MPI processes and will add this to the local + * DynamicSparsityPattern. + * + * @param dsp A dynamic sparsity pattern that has been built locally and which + * we need to extend according to the sparsity of rows stored on other MPI + * processes. + * + * @param owned_rows_per_processor A vector containing owned rows for each + * process in the MPI communicator. This input should be the same on all MPI + * processes. + * + * @param mpi_comm The MPI communicator shared between the processors that + * participate in this operation. + * + * @param ghost_range The range of rows this MPI process needs to gather. + * Only a part which is not included in the locally owned rows will be used. + */ + void + gather_sparsity_pattern(DynamicSparsityPattern & dsp, + const std::vector &owned_rows_per_processor, + const MPI_Comm & mpi_comm, + const IndexSet & ghost_range); + #endif diff --git a/source/lac/dynamic_sparsity_pattern.cc b/source/lac/dynamic_sparsity_pattern.cc index 95ca12cdc4..c697a45a1b 100644 --- a/source/lac/dynamic_sparsity_pattern.cc +++ b/source/lac/dynamic_sparsity_pattern.cc @@ -406,6 +406,25 @@ DynamicSparsityPattern::symmetrize() +void +DynamicSparsityPattern::clear_row(const size_type row) +{ + AssertIndexRange(row, n_rows()); + if (!have_entries) + return; + + if (rowset.size() > 0 && !rowset.is_element(row)) + return; + + const size_type rowindex = + rowset.size() == 0 ? row : rowset.index_within_set(row); + + AssertIndexRange(rowindex, lines.size()); + lines[rowindex].entries = std::vector(); +} + + + template void DynamicSparsityPattern::compute_Tmmult_pattern( diff --git a/source/lac/sparsity_tools.cc b/source/lac/sparsity_tools.cc index 1a18621296..b62a712c26 100644 --- a/source/lac/sparsity_tools.cc +++ b/source/lac/sparsity_tools.cc @@ -911,6 +911,117 @@ namespace SparsityTools #ifdef DEAL_II_WITH_MPI + + void + gather_sparsity_pattern(DynamicSparsityPattern & dsp, + const std::vector &owned, + const MPI_Comm & mpi_comm, + const IndexSet & ghost_range) + { + const unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm); + + AssertDimension(owned.size(), Utilities::MPI::n_mpi_processes(mpi_comm)); +# ifdef DEBUG + for (const auto &set : owned) + Assert(set.is_contiguous(), ExcNotImplemented()); +# endif + + Assert(owned[myid].is_ascending_and_one_to_one(mpi_comm), + ExcNotImplemented()); + + std::vector start_index(owned.size() + + 1); + start_index[0] = 0; + for (DynamicSparsityPattern::size_type i = 0; i < owned.size(); ++i) + start_index[i + 1] = start_index[i] + owned[i].n_elements(); + + using map_vec_t = + std::map>; + + // 1. limit rows to non owned: + IndexSet requested_rows(ghost_range); + requested_rows.subtract_set(owned[myid]); + + // 2. go through requested_rows, figure out the owner and add the row to + // request + map_vec_t rows_data; + { + unsigned int cpu_rank = 0; + for (const auto &row : requested_rows) + { + // calculate destination CPU + while (row >= start_index[cpu_rank + 1]) + ++cpu_rank; + + // since we removed owned, we should not end up with + // our rows + Assert(cpu_rank != myid, ExcInternalError()); + + rows_data[cpu_rank].push_back(row); + } + } + + // 3. get what others ask us to send + const auto rows_data_received = + Utilities::MPI::some_to_some(mpi_comm, rows_data); + + // 4. now prepare data to be sent in the same format as in + // distribute_sparsity_pattern() below, i.e. + // rX,num_rX,cols_rX + map_vec_t send_data; + for (const auto &data : rows_data_received) + { + for (const auto &row : data.second) + { + const auto rlen = dsp.row_length(row); + + // skip empty lines + if (rlen == 0) + continue; + + // save entries + send_data[data.first].push_back(row); // row index + send_data[data.first].push_back(rlen); // number of entries + for (DynamicSparsityPattern::size_type c = 0; c < rlen; ++c) + send_data[data.first].push_back( + dsp.column_number(row, c)); // columns + } // loop over rows + } // loop over received data + + // 5. communicate rows + const auto received_data = + Utilities::MPI::some_to_some(mpi_comm, send_data); + + // 6. add result to our sparsity + for (const auto &data : received_data) + { + const auto &recv_buf = data.second; + auto ptr = recv_buf.begin(); + const auto end = recv_buf.end(); + while (ptr != end) + { + const auto row = *(ptr++); + Assert(ptr != end, ExcInternalError()); + + const auto n_entries = *(ptr++); + Assert(n_entries > 0, ExcInternalError()); + Assert(ptr != end, ExcInternalError()); + + // make sure we clear whatever was previously stored + // in these rows. Otherwise we can't guarantee that the + // data is consistent across MPI communicator. + dsp.clear_row(row); + + Assert(ptr + (n_entries - 1) != end, ExcInternalError()); + dsp.add_entries(row, ptr, ptr + n_entries, true); + ptr += n_entries; + } + Assert(ptr == end, ExcInternalError()); + } + } + + + void distribute_sparsity_pattern( DynamicSparsityPattern & dsp, diff --git a/tests/sparsity/dynamic_sparsity_pattern_19.cc b/tests/sparsity/dynamic_sparsity_pattern_19.cc new file mode 100644 index 0000000000..b0b8d0f497 --- /dev/null +++ b/tests/sparsity/dynamic_sparsity_pattern_19.cc @@ -0,0 +1,105 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check SparsityTools::gather_sparsity_pattern() + +#include + +#include "../tests.h" + + +void +test() +{ + const unsigned int N = 10; + MPI_Comm comm(MPI_COMM_WORLD); + const unsigned int myid = Utilities::MPI::this_mpi_process(comm); + + IndexSet owned(N); + if (myid == 0) + owned.add_range(0, 3); + else if (myid == 1) + owned.add_range(3, 7); + else + owned.add_range(7, 10); + + const auto rows_per_cpu = Utilities::MPI::all_gather(comm, owned); + + IndexSet ghost_range(owned); + if (myid == 0) + { + ghost_range.add_range(3, 5); + } + else if (myid == 1) + { + ghost_range.add_range(1, 3); + ghost_range.add_range(7, 9); + } + else + { + ghost_range.add_range(4, 7); + } + + DynamicSparsityPattern dsp; + dsp.reinit(N, N); + // hard-code some sparsity + if (myid == 0) + { + for (unsigned int i = 0; i < 3; ++i) + for (unsigned int j = 0; j < 5 + i; ++j) + dsp.add(i, j); + + // add to ghost to make sure they are reset + dsp.add(3, 9); + } + else if (myid == 1) + { + for (unsigned int i = 3; i < 5; ++i) + for (unsigned int j = 0; j < 5 + i; ++j) + dsp.add(i, j); + + for (unsigned int i = 5; i < 7; ++i) + for (unsigned int j = i - 4; j < 10; ++j) + dsp.add(i, j); + } + else + { + for (unsigned int i = 7; i < 10; ++i) + for (unsigned int j = i - 4; j < 10; ++j) + dsp.add(i, j); + } + + deallog << "Before gather_sparsity_pattern:" << std::endl; + dsp.print(deallog.get_file_stream()); + + SparsityTools::gather_sparsity_pattern(dsp, rows_per_cpu, comm, ghost_range); + + deallog << "After gather_sparsity_pattern:" << std::endl; + dsp.print(deallog.get_file_stream()); +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + MPILogInitAll log; + + test(); + return 0; +} diff --git a/tests/sparsity/dynamic_sparsity_pattern_19.mpirun=3.output b/tests/sparsity/dynamic_sparsity_pattern_19.mpirun=3.output new file mode 100644 index 0000000000..fd13cba169 --- /dev/null +++ b/tests/sparsity/dynamic_sparsity_pattern_19.mpirun=3.output @@ -0,0 +1,71 @@ + +DEAL:0::Before gather_sparsity_pattern: +[0,0,1,2,3,4] +[1,0,1,2,3,4,5] +[2,0,1,2,3,4,5,6] +[3,9] +[4] +[5] +[6] +[7] +[8] +[9] +DEAL:0::After gather_sparsity_pattern: +[0,0,1,2,3,4] +[1,0,1,2,3,4,5] +[2,0,1,2,3,4,5,6] +[3,0,1,2,3,4,5,6,7] +[4,0,1,2,3,4,5,6,7,8] +[5] +[6] +[7] +[8] +[9] + +DEAL:1::Before gather_sparsity_pattern: +[0] +[1] +[2] +[3,0,1,2,3,4,5,6,7] +[4,0,1,2,3,4,5,6,7,8] +[5,1,2,3,4,5,6,7,8,9] +[6,2,3,4,5,6,7,8,9] +[7] +[8] +[9] +DEAL:1::After gather_sparsity_pattern: +[0] +[1,0,1,2,3,4,5] +[2,0,1,2,3,4,5,6] +[3,0,1,2,3,4,5,6,7] +[4,0,1,2,3,4,5,6,7,8] +[5,1,2,3,4,5,6,7,8,9] +[6,2,3,4,5,6,7,8,9] +[7,3,4,5,6,7,8,9] +[8,4,5,6,7,8,9] +[9] + + +DEAL:2::Before gather_sparsity_pattern: +[0] +[1] +[2] +[3] +[4] +[5] +[6] +[7,3,4,5,6,7,8,9] +[8,4,5,6,7,8,9] +[9,5,6,7,8,9] +DEAL:2::After gather_sparsity_pattern: +[0] +[1] +[2] +[3] +[4,0,1,2,3,4,5,6,7,8] +[5,1,2,3,4,5,6,7,8,9] +[6,2,3,4,5,6,7,8,9] +[7,3,4,5,6,7,8,9] +[8,4,5,6,7,8,9] +[9,5,6,7,8,9] + -- 2.39.5