Reference documentation for deal.II version GIT relicensing-446-ge11b936273 2024-04-20 12:30:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
sparse_matrix_tools.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2022 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_sparse_matrix_tools_h
16#define dealii_sparse_matrix_tools_h
17
18#include <deal.II/base/config.h>
19
21
23
25
27
32{
55 template <typename SparseMatrixType,
56 typename SparsityPatternType,
57 typename SparseMatrixType2,
58 typename SparsityPatternType2>
59 void
60 restrict_to_serial_sparse_matrix(const SparseMatrixType &sparse_matrix_in,
61 const SparsityPatternType &sparsity_pattern,
62 const IndexSet &requested_is,
63 SparseMatrixType2 &system_matrix_out,
64 SparsityPatternType2 &sparsity_pattern_out);
65
78 template <typename SparseMatrixType,
79 typename SparsityPatternType,
80 typename SparseMatrixType2,
81 typename SparsityPatternType2>
82 void
83 restrict_to_serial_sparse_matrix(const SparseMatrixType &sparse_matrix_in,
84 const SparsityPatternType &sparsity_pattern,
85 const IndexSet &index_set_0,
86 const IndexSet &index_set_1,
87 SparseMatrixType2 &system_matrix_out,
88 SparsityPatternType2 &sparsity_pattern_out);
89
105 template <int dim,
106 int spacedim,
107 typename SparseMatrixType,
108 typename SparsityPatternType,
109 typename Number>
110 void
111 restrict_to_cells(const SparseMatrixType &system_matrix,
112 const SparsityPatternType &sparsity_pattern,
113 const DoFHandler<dim, spacedim> &dof_handler,
114 std::vector<FullMatrix<Number>> &blocks);
115
129 template <typename SparseMatrixType,
130 typename SparsityPatternType,
131 typename Number>
132 void
134 const SparseMatrixType &sparse_matrix_in,
135 const SparsityPatternType &sparsity_pattern,
136 const std::vector<std::vector<types::global_dof_index>> &indices,
137 std::vector<FullMatrix<Number>> &blocks);
138
139
140#ifndef DOXYGEN
141 /*---------------------- Inline functions ---------------------------------*/
142
143 namespace internal
144 {
145 template <typename T>
146 using get_mpi_communicator_t =
147 decltype(std::declval<const T>().get_mpi_communicator());
148
149 template <typename T>
150 constexpr bool has_get_mpi_communicator =
151 ::internal::is_supported_operation<get_mpi_communicator_t, T>;
152
153 template <typename T>
154 using local_size_t = decltype(std::declval<const T>().local_size());
155
156 template <typename T>
157 constexpr bool has_local_size =
158 ::internal::is_supported_operation<local_size_t, T>;
159
160 template <typename SparseMatrixType,
161 std::enable_if_t<has_get_mpi_communicator<SparseMatrixType>,
162 SparseMatrixType> * = nullptr>
164 get_mpi_communicator(const SparseMatrixType &sparse_matrix)
165 {
166 return sparse_matrix.get_mpi_communicator();
167 }
168
169 template <typename SparseMatrixType,
170 std::enable_if_t<!has_get_mpi_communicator<SparseMatrixType>,
171 SparseMatrixType> * = nullptr>
173 get_mpi_communicator(const SparseMatrixType &sparse_matrix)
174 {
175 (void)sparse_matrix;
176 return MPI_COMM_SELF;
177 }
178
179 template <typename SparseMatrixType,
180 std::enable_if_t<has_local_size<SparseMatrixType>,
181 SparseMatrixType> * = nullptr>
182 unsigned int
183 get_local_size(const SparseMatrixType &sparse_matrix)
184 {
185 return sparse_matrix.local_size();
186 }
187
188 template <typename SparseMatrixType,
189 std::enable_if_t<!has_local_size<SparseMatrixType>,
190 SparseMatrixType> * = nullptr>
191 unsigned int
192 get_local_size(const SparseMatrixType &sparse_matrix)
193 {
194 AssertDimension(sparse_matrix.m(), sparse_matrix.n());
195
196 return sparse_matrix.m();
197 }
198
199 // Helper function to extract for a distributed sparse matrix rows
200 // potentially not owned by the current process.
201 template <typename Number,
202 typename SparseMatrixType,
203 typename SparsityPatternType>
204 std::vector<std::vector<std::pair<types::global_dof_index, Number>>>
205 extract_remote_rows(const SparseMatrixType &system_matrix,
206 const SparsityPatternType &sparsity_pattern,
207 const IndexSet &locally_active_dofs,
208 const MPI_Comm comm)
209 {
210 std::vector<unsigned int> dummy(locally_active_dofs.n_elements());
211
212 const auto local_size = get_local_size(system_matrix);
213 const auto [prefix_sum, total_sum] =
215 IndexSet locally_owned_dofs(total_sum);
216 locally_owned_dofs.add_range(prefix_sum, prefix_sum + local_size);
217
219 process(locally_owned_dofs, locally_active_dofs, comm, dummy, true);
220
222 std::vector<
223 std::pair<types::global_dof_index, types::global_dof_index>>,
224 std::vector<unsigned int>>
225 consensus_algorithm;
226 consensus_algorithm.run(process, comm);
227
228 using T1 = std::vector<
229 std::pair<types::global_dof_index,
230 std::vector<std::pair<types::global_dof_index, Number>>>>;
231
232 auto requesters = process.get_requesters();
233
234 std::vector<std::vector<std::pair<types::global_dof_index, Number>>>
235 locally_relevant_matrix_entries(locally_active_dofs.n_elements());
236
237
238 std::vector<unsigned int> ranks;
239 ranks.reserve(requesters.size());
240
241 for (const auto &i : requesters)
242 ranks.push_back(i.first);
243
244 std::vector<std::vector<unsigned int>> row_to_procs(
245 locally_owned_dofs.n_elements());
246
247 for (const auto &requester : requesters)
248 for (const auto &index : requester.second)
249 row_to_procs[locally_owned_dofs.index_within_set(index)].push_back(
250 requester.first);
251
252 std::map<unsigned int, T1> data;
253
254 std::pair<types::global_dof_index,
255 std::vector<std::pair<types::global_dof_index, Number>>>
256 buffer;
257
258 for (unsigned int i = 0; i < row_to_procs.size(); ++i)
259 {
260 if (row_to_procs[i].empty())
261 continue;
262
263 const auto row = locally_owned_dofs.nth_index_in_set(i);
264 auto entry = system_matrix.begin(row);
265
266 const unsigned int row_length = sparsity_pattern.row_length(row);
267
268 buffer.first = row;
269 buffer.second.resize(row_length);
270
271 for (unsigned int j = 0; j < row_length; ++j, ++entry)
272 buffer.second[j] = {entry->column(), entry->value()};
273
274 for (const auto &proc : row_to_procs[i])
275 data[proc].emplace_back(buffer);
276 }
277
278 ::Utilities::MPI::ConsensusAlgorithms::selector<T1>(
279 ranks,
280 [&](const unsigned int other_rank) { return data[other_rank]; },
281 [&](const unsigned int &, const T1 &buffer_recv) {
282 for (const auto &i : buffer_recv)
283 {
284 auto &dst =
285 locally_relevant_matrix_entries[locally_active_dofs
286 .index_within_set(i.first)];
287 dst = i.second;
288 std::sort(dst.begin(),
289 dst.end(),
290 [](const auto &a, const auto &b) {
291 return a.first < b.first;
292 });
293 }
294 },
295 comm);
296
297 return locally_relevant_matrix_entries;
298 }
299 } // namespace internal
300
301
302
303 template <typename SparseMatrixType,
304 typename SparsityPatternType,
305 typename SparseMatrixType2,
306 typename SparsityPatternType2>
307 void
308 restrict_to_serial_sparse_matrix(const SparseMatrixType &system_matrix,
309 const SparsityPatternType &sparsity_pattern,
310 const IndexSet &index_set_0,
311 const IndexSet &index_set_1,
312 SparseMatrixType2 &system_matrix_out,
313 SparsityPatternType2 &sparsity_pattern_out)
314 {
315 Assert(index_set_1.size() == 0 || index_set_0.size() == index_set_1.size(),
317
318 auto index_set_1_cleared = index_set_1;
319 if (index_set_1.size() != 0)
320 index_set_1_cleared.subtract_set(index_set_0);
321
322 const auto index_within_set = [&index_set_0,
323 &index_set_1_cleared](const auto n) {
324 if (index_set_0.is_element(n))
325 return index_set_0.index_within_set(n);
326 else
327 return index_set_0.n_elements() +
328 index_set_1_cleared.index_within_set(n);
329 };
330
331 // 1) collect needed rows
332 auto index_set_union = index_set_0;
333
334 if (index_set_1.size() != 0)
335 index_set_union.add_indices(index_set_1_cleared);
336
337 // TODO: actually only communicate remote rows as in the case of
338 // SparseMatrixTools::restrict_to_cells()
339 const auto locally_relevant_matrix_entries =
340 internal::extract_remote_rows<typename SparseMatrixType2::value_type>(
341 system_matrix,
342 sparsity_pattern,
343 index_set_union,
344 internal::get_mpi_communicator(system_matrix));
345
346
347 // 2) create sparsity pattern
348 const unsigned int n_rows = index_set_union.n_elements();
349 const unsigned int n_cols = index_set_union.n_elements();
350 const unsigned int entries_per_row =
351 locally_relevant_matrix_entries.empty() ?
352 0 :
353 std::max_element(locally_relevant_matrix_entries.begin(),
354 locally_relevant_matrix_entries.end(),
355 [](const auto &a, const auto &b) {
356 return a.size() < b.size();
357 })
358 ->size();
359
360 sparsity_pattern_out.reinit(n_rows, n_cols, entries_per_row);
361
362 std::vector<types::global_dof_index> temp_indices;
363 std::vector<typename SparseMatrixType2::value_type> temp_values;
364
365 for (unsigned int row = 0; row < index_set_union.n_elements(); ++row)
366 {
367 const auto &global_row_entries = locally_relevant_matrix_entries[row];
368
369 temp_indices.clear();
370 temp_indices.reserve(global_row_entries.size());
371
372 for (const auto &global_row_entry : global_row_entries)
373 {
374 const auto global_index = std::get<0>(global_row_entry);
375
376 if (index_set_union.is_element(global_index))
377 temp_indices.push_back(index_within_set(global_index));
378 }
379
380 sparsity_pattern_out.add_entries(
381 index_within_set(index_set_union.nth_index_in_set(row)),
382 temp_indices.begin(),
383 temp_indices.end());
384 }
385
386 sparsity_pattern_out.compress();
387
388 // 3) setup matrix
389 system_matrix_out.reinit(sparsity_pattern_out);
390
391 // 4) fill entries
392 for (unsigned int row = 0; row < index_set_union.n_elements(); ++row)
393 {
394 const auto &global_row_entries = locally_relevant_matrix_entries[row];
395
396 temp_indices.clear();
397 temp_values.clear();
398
399 temp_indices.reserve(global_row_entries.size());
400 temp_values.reserve(global_row_entries.size());
401
402 for (const auto &global_row_entry : global_row_entries)
403 {
404 const auto global_index = std::get<0>(global_row_entry);
405
406 if (index_set_union.is_element(global_index))
407 {
408 temp_indices.push_back(index_within_set(global_index));
409 temp_values.push_back(std::get<1>(global_row_entry));
410 }
411 }
412
413 system_matrix_out.add(index_within_set(
414 index_set_union.nth_index_in_set(row)),
415 temp_indices,
416 temp_values);
417 }
418
419 system_matrix_out.compress(VectorOperation::add);
420 }
421
422
423
424 template <typename SparseMatrixType,
425 typename SparsityPatternType,
426 typename SparseMatrixType2,
427 typename SparsityPatternType2>
428 void
429 restrict_to_serial_sparse_matrix(const SparseMatrixType &system_matrix,
430 const SparsityPatternType &sparsity_pattern,
431 const IndexSet &requested_is,
432 SparseMatrixType2 &system_matrix_out,
433 SparsityPatternType2 &sparsity_pattern_out)
434 {
436 sparsity_pattern,
437 requested_is,
438 IndexSet(), // simply pass empty index set
439 system_matrix_out,
440 sparsity_pattern_out);
441 }
442
443
444
445 template <typename SparseMatrixType,
446 typename SparsityPatternType,
447 typename Number>
448 void
450 const SparseMatrixType &system_matrix,
451 const SparsityPatternType &sparsity_pattern,
452 const std::vector<std::vector<types::global_dof_index>> &indices,
453 std::vector<FullMatrix<Number>> &blocks)
454 {
455 // 0) determine which rows are locally owned and which ones are remote
456 const auto local_size = internal::get_local_size(system_matrix);
457 const auto prefix_sum = Utilities::MPI::partial_and_total_sum(
458 local_size, internal::get_mpi_communicator(system_matrix));
459 IndexSet locally_owned_dofs(std::get<1>(prefix_sum));
460 locally_owned_dofs.add_range(std::get<0>(prefix_sum),
461 std::get<0>(prefix_sum) + local_size);
462
463 std::vector<::types::global_dof_index> ghost_indices_vector;
464
465 for (const auto &i : indices)
466 ghost_indices_vector.insert(ghost_indices_vector.end(),
467 i.begin(),
468 i.end());
469
470 std::sort(ghost_indices_vector.begin(), ghost_indices_vector.end());
471
472 IndexSet locally_active_dofs(std::get<1>(prefix_sum));
473 locally_active_dofs.add_indices(ghost_indices_vector.begin(),
474 ghost_indices_vector.end());
475
476 locally_active_dofs.subtract_set(locally_owned_dofs);
477
478 // 1) collect remote rows of sparse matrix
479 const auto locally_relevant_matrix_entries =
480 internal::extract_remote_rows<Number>(system_matrix,
481 sparsity_pattern,
482 locally_active_dofs,
483 internal::get_mpi_communicator(
484 system_matrix));
485
486
487 // 2) loop over all cells and "revert" assembly
488 blocks.clear();
489 blocks.resize(indices.size());
490
491 for (unsigned int c = 0; c < indices.size(); ++c)
492 {
493 if (indices[c].empty())
494 continue;
495
496 const auto &local_dof_indices = indices[c];
497 auto &cell_matrix = blocks[c];
498
499 // allocate memory
500 const unsigned int dofs_per_cell = indices[c].size();
501
502 cell_matrix = FullMatrix<Number>(dofs_per_cell, dofs_per_cell);
503
504 // loop over all entries of the restricted element matrix and
505 // do different things if rows are locally owned or not and
506 // if column entries of that row exist or not
507 for (unsigned int i = 0; i < dofs_per_cell; ++i)
508 for (unsigned int j = 0; j < dofs_per_cell; ++j)
509 {
510 if (locally_owned_dofs.is_element(
511 local_dof_indices[i])) // row is local
512 {
513 cell_matrix(i, j) =
514 sparsity_pattern.exists(local_dof_indices[i],
515 local_dof_indices[j]) ?
516 system_matrix(local_dof_indices[i],
517 local_dof_indices[j]) :
518 0;
519 }
520 else // row is ghost
521 {
522 Assert(locally_active_dofs.is_element(local_dof_indices[i]),
524
525 const auto &row_entries =
526 locally_relevant_matrix_entries[locally_active_dofs
528 local_dof_indices[i])];
529
530 const auto ptr =
531 std::lower_bound(row_entries.begin(),
532 row_entries.end(),
533 std::pair<types::global_dof_index, Number>{
534 local_dof_indices[j], /*dummy*/ 0.0},
535 [](const auto a, const auto b) {
536 return a.first < b.first;
537 });
538
539 if (ptr != row_entries.end() &&
540 local_dof_indices[j] == ptr->first)
541 cell_matrix(i, j) = ptr->second;
542 else
543 cell_matrix(i, j) = 0.0;
544 }
545 }
546 }
547 }
548
549
550
551 template <int dim,
552 int spacedim,
553 typename SparseMatrixType,
554 typename SparsityPatternType,
555 typename Number>
556 void
557 restrict_to_cells(const SparseMatrixType &system_matrix,
558 const SparsityPatternType &sparsity_pattern,
559 const DoFHandler<dim, spacedim> &dof_handler,
560 std::vector<FullMatrix<Number>> &blocks)
561 {
562 std::vector<std::vector<types::global_dof_index>> all_dof_indices;
563 all_dof_indices.resize(dof_handler.get_triangulation().n_active_cells());
564
565 for (const auto &cell : dof_handler.active_cell_iterators())
566 {
567 if (cell->is_locally_owned() == false)
568 continue;
569
570 auto &local_dof_indices = all_dof_indices[cell->active_cell_index()];
571 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
572 cell->get_dof_indices(local_dof_indices);
573 }
574
575 restrict_to_full_matrices(system_matrix,
576 sparsity_pattern,
577 all_dof_indices,
578 blocks);
579 }
580#endif
581
582} // namespace SparseMatrixTools
583
585
586#endif
const Triangulation< dim, spacedim > & get_triangulation() const
size_type size() const
Definition index_set.h:1765
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1990
size_type n_elements() const
Definition index_set.h:1923
bool is_element(const size_type index) const
Definition index_set.h:1883
void subtract_set(const IndexSet &other)
Definition index_set.cc:470
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1820
unsigned int n_active_cells() const
virtual std::vector< unsigned int > run(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm) override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
void restrict_to_cells(const SparseMatrixType &system_matrix, const SparsityPatternType &sparsity_pattern, const DoFHandler< dim, spacedim > &dof_handler, std::vector< FullMatrix< Number > > &blocks)
void restrict_to_serial_sparse_matrix(const SparseMatrixType &sparse_matrix_in, const SparsityPatternType &sparsity_pattern, const IndexSet &requested_is, SparseMatrixType2 &system_matrix_out, SparsityPatternType2 &sparsity_pattern_out)
void restrict_to_full_matrices(const SparseMatrixType &sparse_matrix_in, const SparsityPatternType &sparsity_pattern, const std::vector< std::vector< types::global_dof_index > > &indices, std::vector< FullMatrix< Number > > &blocks)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
std::pair< T, T > partial_and_total_sum(const T &value, const MPI_Comm comm)
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const Iterator & begin
Definition parallel.h:609
unsigned int global_dof_index
Definition types.h:81
*braid_SplitCommworld & comm