--- /dev/null
+New: Add SparsityTools::gather_sparsity_pattern().
+<br>
+(Denis Davydov, 2019/03/19)
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.
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<IndexSet> &owned_rows_per_processor,
+ const MPI_Comm & mpi_comm,
+ const IndexSet & ghost_range);
+
#endif
+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<size_type>();
+}
+
+
+
template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
void
DynamicSparsityPattern::compute_Tmmult_pattern(
#ifdef DEAL_II_WITH_MPI
+
+ void
+ gather_sparsity_pattern(DynamicSparsityPattern & dsp,
+ const std::vector<IndexSet> &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<DynamicSparsityPattern::size_type> 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<unsigned int, std::vector<DynamicSparsityPattern::size_type>>;
+
+ // 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,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/sparsity_tools.h>
+
+#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;
+}
--- /dev/null
+
+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]
+
projects / dealii.git / commitdiff