10#include <deal.II/base/vectorization.h>
11#include <deal.II/lac/sparse_matrix.h>
16 template <
int simd_length>
20 , mpi_communicator(MPI_COMM_SELF)
25 template <
int simd_length>
27 const unsigned int n_internal_dofs,
28 const dealii::DynamicSparsityPattern &sparsity,
29 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
32 , mpi_communicator(MPI_COMM_SELF)
38 template <
int simd_length>
40 const unsigned int n_internal_dofs,
41 const dealii::DynamicSparsityPattern &dsp,
42 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
45 this->mpi_communicator = partitioner->get_mpi_communicator();
47 this->n_internal_dofs = n_internal_dofs;
48 this->n_locally_owned_dofs = partitioner->locally_owned_size();
49 this->partitioner = partitioner;
51 const auto n_locally_relevant_dofs =
52 partitioner->locally_owned_size() + partitioner->n_ghost_indices();
61 dealii::DynamicSparsityPattern dsp_minimal(n_locally_relevant_dofs,
62 n_locally_relevant_dofs);
63 for (
unsigned int i = 0; i < n_locally_owned_dofs; ++i) {
64 const auto global_row = partitioner->local_to_global(i);
65 for (
auto it = dsp.begin(global_row); it != dsp.end(global_row); ++it) {
66 const auto global_column = it->column();
67 const auto j = partitioner->global_to_local(global_column);
68 dsp_minimal.add(i, j);
69 if (j >= n_locally_owned_dofs) {
70 Assert(j < n_locally_relevant_dofs, dealii::ExcInternalError());
71 dsp_minimal.add(j, i);
76 dealii::SparsityPattern sparsity;
77 sparsity.copy_from(dsp_minimal);
79 Assert(n_internal_dofs <= sparsity.n_rows(), dealii::ExcInternalError());
80 Assert(n_internal_dofs % simd_length == 0, dealii::ExcInternalError());
81 Assert(n_internal_dofs <= n_locally_owned_dofs, dealii::ExcInternalError());
82 Assert(n_locally_owned_dofs <= sparsity.n_rows(),
83 dealii::ExcInternalError());
85 row_starts.resize_fast(sparsity.n_rows() + 1);
86 column_indices.resize_fast(sparsity.n_nonzero_elements());
87 indices_transposed.resize_fast(sparsity.n_nonzero_elements());
88 AssertThrow(sparsity.n_nonzero_elements() <
89 std::numeric_limits<unsigned int>::max(),
90 dealii::ExcMessage(
"Transposed indices only support up to 4 "
91 "billion matrix entries per MPI rank. Try to"
92 " split into smaller problems with MPI"));
98 unsigned int *col_ptr = column_indices.data();
99 unsigned int *transposed_ptr = indices_transposed.data();
101 for (
unsigned int i = 0; i < n_internal_dofs; i += simd_length) {
102 auto jts = generate_iterators<simd_length>(
103 [&](
auto k) {
return sparsity.begin(i + k); });
106 for (
unsigned int k = 0; k < simd_length; ++k) {
107 const unsigned int column = jts[k]->column();
109 const std::size_t position = sparsity(column, i + k);
110 if (column < n_internal_dofs) {
111 const unsigned int my_row_length = sparsity.row_length(column);
112 const std::size_t position_diag = sparsity(column, column);
113 const std::size_t pos_within_row = position - position_diag;
114 const unsigned int simd_offset = column % simd_length;
115 *transposed_ptr++ = position - simd_offset * my_row_length -
116 pos_within_row + simd_offset +
117 pos_within_row * simd_length;
119 *transposed_ptr++ = position;
122 row_starts[i / simd_length + 1] = col_ptr - column_indices.data();
127 row_starts[n_internal_dofs] = row_starts[n_internal_dofs / simd_length];
129 for (
unsigned int i = n_internal_dofs; i < sparsity.n_rows(); ++i) {
130 for (
auto j = sparsity.begin(i); j != sparsity.end(i); ++j) {
131 const unsigned int column = j->column();
133 const std::size_t position = sparsity(column, i);
134 if (column < n_internal_dofs) {
135 const unsigned int my_row_length = sparsity.row_length(column);
136 const std::size_t position_diag = sparsity(column, column);
137 const std::size_t pos_within_row = position - position_diag;
138 const unsigned int simd_offset = column % simd_length;
139 *transposed_ptr++ = position - simd_offset * my_row_length -
140 pos_within_row + simd_offset +
141 pos_within_row * simd_length;
143 *transposed_ptr++ = position;
145 row_starts[i + 1] = col_ptr - column_indices.data();
148 Assert(col_ptr == column_indices.end(), dealii::ExcInternalError());
152 if (sparsity.n_rows() > n_locally_owned_dofs) {
164 const auto &ghost_targets = partitioner->ghost_targets();
166 receive_targets.resize(ghost_targets.size());
168 for (
unsigned int p = 0; p < receive_targets.size(); ++p) {
169 receive_targets[p].first = ghost_targets[p].first;
172 const auto gt_begin = ghost_targets.begin();
173 auto gt_ptr = ghost_targets.begin();
174 std::size_t index = 0;
175 unsigned int row_count = 0;
177 for (
unsigned int i = n_locally_owned_dofs; i < sparsity.n_rows(); ++i) {
178 index += sparsity.row_length(i);
180 if (row_count == gt_ptr->second) {
181 receive_targets[gt_ptr - gt_begin].second = index;
187 Assert(gt_ptr == partitioner->ghost_targets().end(),
188 dealii::ExcInternalError());
200 std::vector<unsigned int> ghost_ranges(ghost_targets.size() + 1);
201 ghost_ranges[0] = n_locally_owned_dofs;
202 for (
unsigned int p = 0; p < receive_targets.size(); ++p) {
203 ghost_ranges[p + 1] = ghost_ranges[p] + ghost_targets[p].second;
206 std::vector<unsigned int> import_indices_part;
207 for (
auto i : partitioner->import_indices())
208 for (
unsigned int j = i.first; j < i.second; ++j)
209 import_indices_part.push_back(j);
211 AssertDimension(import_indices_part.size(),
212 partitioner->n_import_indices());
214 const auto &import_targets = partitioner->import_targets();
215 entries_to_be_sent.clear();
216 send_targets.resize(import_targets.size());
217 auto idx = import_indices_part.begin();
224 unsigned int p_match = 0;
225 for (
unsigned int p = 0; p < import_targets.size(); ++p) {
237 p_match = (p_match == receive_targets.size() ? 0 : p_match);
238 while (p_match < receive_targets.size() &&
239 receive_targets[p_match].first != import_targets[p].first)
242 for (
unsigned int c = 0; c < import_targets[p].second; ++c, ++idx) {
248 if (p_match == receive_targets.size())
251 const unsigned int row = *idx;
253 entries_to_be_sent.emplace_back(row, 0);
254 for (
auto jt = ++sparsity.begin(row); jt != sparsity.end(row); ++jt) {
255 if (jt->column() >= ghost_ranges[p_match] &&
256 jt->column() < ghost_ranges[p_match + 1]) {
257 const unsigned int position_within_column =
258 jt - sparsity.begin(row);
259 entries_to_be_sent.emplace_back(row, position_within_column);
264 send_targets[p].first = partitioner->import_targets()[p].first;
265 send_targets[p].second = entries_to_be_sent.size();
271 template <
typename Number,
int n_components,
int simd_length>
278 template <
typename Number,
int n_components,
int simd_length>
281 : sparsity(&sparsity)
283 data.resize(
sparsity.n_nonzero_elements() * n_components);
287 template <
typename Number,
int n_components,
int simd_length>
291 this->sparsity = &sparsity;
296 template <
typename Number,
int n_components,
int simd_length>
297 template <
typename SparseMatrix>
299 const std::array<SparseMatrix, n_components> &sparse_matrix,
300 bool locally_indexed )
312 for (
unsigned int i = 0; i < sparsity->n_internal_dofs; i += simd_length) {
314 const unsigned int row_length = sparsity->row_length(i);
316 const unsigned int *js = sparsity->columns(i);
317 for (
unsigned int col_idx = 0; col_idx < row_length;
318 ++col_idx, js += simd_length) {
320 dealii::Tensor<1, n_components, VectorizedArray> temp;
321 for (
unsigned int k = 0; k < simd_length; ++k)
322 for (
unsigned int d = 0; d < n_components; ++d)
324 temp[d][k] = sparse_matrix[d](i + k, js[k]);
326 temp[d][k] = sparse_matrix[d].el(
327 sparsity->partitioner->local_to_global(i + k),
328 sparsity->partitioner->local_to_global(js[k]));
335 for (
unsigned int i = sparsity->n_internal_dofs;
336 i < sparsity->n_locally_owned_dofs;
338 const unsigned int row_length = sparsity->row_length(i);
339 const unsigned int *js = sparsity->columns(i);
340 for (
unsigned int col_idx = 0; col_idx < row_length; ++col_idx, ++js) {
342 dealii::Tensor<1, n_components, Number> temp;
343 for (
unsigned int d = 0; d < n_components; ++d)
345 temp[d] = sparse_matrix[d](i, js[0]);
347 temp[d] = sparse_matrix[d].el(
348 sparsity->partitioner->local_to_global(i),
349 sparsity->partitioner->local_to_global(js[0]));
358 template <
typename Number,
int n_components,
int simd_length>
359 template <
typename SparseMatrix>
361 const SparseMatrix &sparse_matrix,
bool locally_indexed )
373 for (
unsigned int i = 0; i < sparsity->n_internal_dofs; i += simd_length) {
375 const unsigned int row_length = sparsity->row_length(i);
377 const unsigned int *js = sparsity->columns(i);
378 for (
unsigned int col_idx = 0; col_idx < row_length;
379 ++col_idx, js += simd_length) {
382 for (
unsigned int k = 0; k < simd_length; ++k)
384 temp[k] = sparse_matrix(i + k, js[k]);
387 sparse_matrix.el(sparsity->partitioner->local_to_global(i + k),
388 sparsity->partitioner->local_to_global(js[k]));
395 for (
unsigned int i = sparsity->n_internal_dofs;
396 i < sparsity->n_locally_owned_dofs;
399 const unsigned int row_length = sparsity->row_length(i);
400 const unsigned int *js = sparsity->columns(i);
401 for (
unsigned int col_idx = 0; col_idx < row_length; ++col_idx, ++js) {
405 ? sparse_matrix(i, js[0])
407 sparsity->partitioner->local_to_global(i),
408 sparsity->partitioner->local_to_global(js[0]));
dealii::VectorizedArray< Number, simd_length > VectorizedArray
void read_in(const std::array< SparseMatrix, n_components > &sparse_matrix, bool locally_indexed=true)
void reinit(const SparsityPatternSIMD< simd_length > &sparsity)
dealii::AlignedVector< Number > data
const SparsityPatternSIMD< simd_length > * sparsity
void reinit(const unsigned int n_internal_dofs, const dealii::DynamicSparsityPattern &sparsity, const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &partitioner)
unsigned int n_internal_dofs
std::shared_ptr< const dealii::Utilities::MPI::Partitioner > partitioner
std::size_t n_nonzero_elements() const
#define RYUJIN_PARALLEL_REGION_BEGIN
#define RYUJIN_PARALLEL_REGION_END
DEAL_II_ALWAYS_INLINE void increment_iterators(T &iterators)
DEAL_II_ALWAYS_INLINE void write_entry(V &vector, const T &values, unsigned int i)