Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2022 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_sparse_matrix_tools_h
17#define dealii_sparse_matrix_tools_h
18
19#include <deal.II/base/config.h>
20
22
24
26
28
33{
56 template <typename SparseMatrixType,
57 typename SparsityPatternType,
58 typename SparseMatrixType2,
59 typename SparsityPatternType2>
60 void
61 restrict_to_serial_sparse_matrix(const SparseMatrixType & sparse_matrix_in,
62 const SparsityPatternType &sparsity_pattern,
63 const IndexSet & requested_is,
64 SparseMatrixType2 & system_matrix_out,
65 SparsityPatternType2 &sparsity_pattern_out);
66
79 template <typename SparseMatrixType,
80 typename SparsityPatternType,
81 typename SparseMatrixType2,
82 typename SparsityPatternType2>
83 void
84 restrict_to_serial_sparse_matrix(const SparseMatrixType & sparse_matrix_in,
85 const SparsityPatternType &sparsity_pattern,
86 const IndexSet & index_set_0,
87 const IndexSet & index_set_1,
88 SparseMatrixType2 & system_matrix_out,
89 SparsityPatternType2 &sparsity_pattern_out);
90
106 template <int dim,
107 int spacedim,
108 typename SparseMatrixType,
109 typename SparsityPatternType,
110 typename Number>
111 void
112 restrict_to_cells(const SparseMatrixType & system_matrix,
113 const SparsityPatternType & sparsity_pattern,
114 const DoFHandler<dim, spacedim> &dof_handler,
115 std::vector<FullMatrix<Number>> &blocks);
116
130 template <typename SparseMatrixType,
131 typename SparsityPatternType,
132 typename Number>
133 void
135 const SparseMatrixType & sparse_matrix_in,
136 const SparsityPatternType & sparsity_pattern,
137 const std::vector<std::vector<types::global_dof_index>> &indices,
138 std::vector<FullMatrix<Number>> & blocks);
139
140
141#ifndef DOXYGEN
142 /*---------------------- Inline functions ---------------------------------*/
143
144 namespace internal
145 {
146 template <typename T>
147 std::tuple<T, T>
148 compute_prefix_sum(const T &value, const MPI_Comm comm)
149 {
150# ifndef DEAL_II_WITH_MPI
151 (void)comm;
152 return {0, value};
153# else
154 T prefix = {};
155
156 int ierr =
157 MPI_Exscan(&value,
158 &prefix,
159 1,
161 MPI_SUM,
162 comm);
163 AssertThrowMPI(ierr);
164
165 T sum = Utilities::MPI::sum(value, comm);
166
167 return {prefix, sum};
168# endif
169 }
170
171 template <typename T>
172 using get_mpi_communicator_t =
173 decltype(std::declval<T const>().get_mpi_communicator());
174
175 template <typename T>
176 constexpr bool has_get_mpi_communicator =
177 ::internal::is_supported_operation<get_mpi_communicator_t, T>;
178
179 template <typename T>
180 using local_size_t = decltype(std::declval<T const>().local_size());
181
182 template <typename T>
183 constexpr bool has_local_size =
184 ::internal::is_supported_operation<local_size_t, T>;
185
186 template <typename SparseMatrixType,
187 std::enable_if_t<has_get_mpi_communicator<SparseMatrixType>,
188 SparseMatrixType> * = nullptr>
190 get_mpi_communicator(const SparseMatrixType &sparse_matrix)
191 {
192 return sparse_matrix.get_mpi_communicator();
193 }
194
195 template <typename SparseMatrixType,
196 std::enable_if_t<!has_get_mpi_communicator<SparseMatrixType>,
197 SparseMatrixType> * = nullptr>
199 get_mpi_communicator(const SparseMatrixType &sparse_matrix)
200 {
201 (void)sparse_matrix;
202 return MPI_COMM_SELF;
203 }
204
205 template <typename SparseMatrixType,
206 std::enable_if_t<has_local_size<SparseMatrixType>,
207 SparseMatrixType> * = nullptr>
208 unsigned int
209 get_local_size(const SparseMatrixType &sparse_matrix)
210 {
211 return sparse_matrix.local_size();
212 }
213
214 template <typename SparseMatrixType,
215 std::enable_if_t<!has_local_size<SparseMatrixType>,
216 SparseMatrixType> * = nullptr>
217 unsigned int
218 get_local_size(const SparseMatrixType &sparse_matrix)
219 {
220 AssertDimension(sparse_matrix.m(), sparse_matrix.n());
221
222 return sparse_matrix.m();
223 }
224
225 // Helper function to extract for a distributed sparse matrix rows
226 // potentially not owned by the current process.
227 template <typename Number,
228 typename SparseMatrixType,
229 typename SparsityPatternType>
230 std::vector<std::vector<std::pair<types::global_dof_index, Number>>>
231 extract_remote_rows(const SparseMatrixType & system_matrix,
232 const SparsityPatternType &sparsity_pattern,
233 const IndexSet & locally_active_dofs,
234 const MPI_Comm comm)
235 {
236 std::vector<unsigned int> dummy(locally_active_dofs.n_elements());
237
238 const auto local_size = get_local_size(system_matrix);
239 const auto prefix_sum = compute_prefix_sum(local_size, comm);
240 IndexSet locally_owned_dofs(std::get<1>(prefix_sum));
241 locally_owned_dofs.add_range(std::get<0>(prefix_sum),
242 std::get<0>(prefix_sum) + local_size);
243
245 process(locally_owned_dofs, locally_active_dofs, comm, dummy, true);
246
248 std::vector<
249 std::pair<types::global_dof_index, types::global_dof_index>>,
250 std::vector<unsigned int>>
251 consensus_algorithm;
252 consensus_algorithm.run(process, comm);
253
254 using T1 = std::vector<
255 std::pair<types::global_dof_index,
256 std::vector<std::pair<types::global_dof_index, Number>>>>;
257
258 auto requesters = process.get_requesters();
259
260 std::vector<std::vector<std::pair<types::global_dof_index, Number>>>
261 locally_relevant_matrix_entries(locally_active_dofs.n_elements());
262
263
264 std::vector<unsigned int> ranks;
265 ranks.reserve(requesters.size());
266
267 for (const auto &i : requesters)
268 ranks.push_back(i.first);
269
270 std::vector<std::vector<unsigned int>> row_to_procs(
271 locally_owned_dofs.n_elements());
272
273 for (const auto &requester : requesters)
274 for (const auto &index : requester.second)
275 row_to_procs[locally_owned_dofs.index_within_set(index)].push_back(
276 requester.first);
277
278 std::map<unsigned int, T1> data;
279
280 for (unsigned int i = 0; i < row_to_procs.size(); ++i)
281 {
282 if (row_to_procs[i].size() == 0)
283 continue;
284
285 const auto row = locally_owned_dofs.nth_index_in_set(i);
286 auto entry = system_matrix.begin(row);
287
288 const unsigned int row_length = sparsity_pattern.row_length(row);
289
290
291 std::pair<types::global_dof_index,
292 std::vector<std::pair<types::global_dof_index, Number>>>
293 buffer;
294 buffer.first = row;
295
296 for (unsigned int i = 0; i < row_length; ++i)
297 {
298 buffer.second.emplace_back(entry->column(), entry->value());
299
300 if (i + 1 != row_length)
301 ++entry;
302 }
303
304 for (const auto &proc :
305 row_to_procs[locally_owned_dofs.index_within_set(buffer.first)])
306 data[proc].emplace_back(buffer);
307 }
308
309 ::Utilities::MPI::ConsensusAlgorithms::selector<T1>(
310 ranks,
311 [&](const unsigned int other_rank) { return data[other_rank]; },
312 [&](const unsigned int &, const T1 &buffer_recv) {
313 for (const auto &i : buffer_recv)
314 {
315 auto &dst =
316 locally_relevant_matrix_entries[locally_active_dofs
317 .index_within_set(i.first)];
318 dst = i.second;
319 std::sort(dst.begin(),
320 dst.end(),
321 [](const auto &a, const auto &b) {
322 return a.first < b.first;
323 });
324 }
325 },
326 comm);
327
328 return locally_relevant_matrix_entries;
329 }
330 } // namespace internal
331
332
333
334 template <typename SparseMatrixType,
335 typename SparsityPatternType,
336 typename SparseMatrixType2,
337 typename SparsityPatternType2>
338 void
339 restrict_to_serial_sparse_matrix(const SparseMatrixType & system_matrix,
340 const SparsityPatternType &sparsity_pattern,
341 const IndexSet & index_set_0,
342 const IndexSet & index_set_1,
343 SparseMatrixType2 & system_matrix_out,
344 SparsityPatternType2 &sparsity_pattern_out)
345 {
346 Assert(index_set_1.size() == 0 || index_set_0.size() == index_set_1.size(),
348
349 auto index_set_1_cleared = index_set_1;
350 if (index_set_1.size() != 0)
351 index_set_1_cleared.subtract_set(index_set_0);
352
353 const auto index_within_set = [&index_set_0,
354 &index_set_1_cleared](const auto n) {
355 if (index_set_0.is_element(n))
356 return index_set_0.index_within_set(n);
357 else
358 return index_set_0.n_elements() +
359 index_set_1_cleared.index_within_set(n);
360 };
361
362 // 1) collect needed rows
363 auto index_set_union = index_set_0;
364
365 if (index_set_1.size() != 0)
366 index_set_union.add_indices(index_set_1_cleared);
367
368 // TODO: actually only communicate remote rows as in the case of
369 // SparseMatrixTools::restrict_to_cells()
370 const auto locally_relevant_matrix_entries =
371 internal::extract_remote_rows<typename SparseMatrixType2::value_type>(
372 system_matrix,
373 sparsity_pattern,
374 index_set_union,
375 internal::get_mpi_communicator(system_matrix));
376
377
378 // 2) create sparsity pattern
379 const unsigned int n_rows = index_set_union.n_elements();
380 const unsigned int n_cols = index_set_union.n_elements();
381 const unsigned int entries_per_row =
382 locally_relevant_matrix_entries.size() == 0 ?
383 0 :
384 std::max_element(locally_relevant_matrix_entries.begin(),
385 locally_relevant_matrix_entries.end(),
386 [](const auto &a, const auto &b) {
387 return a.size() < b.size();
388 })
389 ->size();
390
391 sparsity_pattern_out.reinit(n_rows, n_cols, entries_per_row);
392
393 std::vector<types::global_dof_index> temp_indices;
394 std::vector<typename SparseMatrixType2::value_type> temp_values;
395
396 for (unsigned int row = 0; row < index_set_union.n_elements(); ++row)
397 {
398 const auto &global_row_entries = locally_relevant_matrix_entries[row];
399
400 temp_indices.clear();
401 temp_indices.reserve(global_row_entries.size());
402
403 for (const auto &global_row_entry : global_row_entries)
404 {
405 const auto global_index = std::get<0>(global_row_entry);
406
407 if (index_set_union.is_element(global_index))
408 temp_indices.push_back(index_within_set(global_index));
409 }
410
411 sparsity_pattern_out.add_entries(
412 index_within_set(index_set_union.nth_index_in_set(row)),
413 temp_indices.begin(),
414 temp_indices.end());
415 }
416
417 sparsity_pattern_out.compress();
418
419 // 3) setup matrix
420 system_matrix_out.reinit(sparsity_pattern_out);
421
422 // 4) fill entries
423 for (unsigned int row = 0; row < index_set_union.n_elements(); ++row)
424 {
425 const auto &global_row_entries = locally_relevant_matrix_entries[row];
426
427 temp_indices.clear();
428 temp_values.clear();
429
430 temp_indices.reserve(global_row_entries.size());
431 temp_values.reserve(global_row_entries.size());
432
433 for (const auto &global_row_entry : global_row_entries)
434 {
435 const auto global_index = std::get<0>(global_row_entry);
436
437 if (index_set_union.is_element(global_index))
438 {
439 temp_indices.push_back(index_within_set(global_index));
440 temp_values.push_back(std::get<1>(global_row_entry));
441 }
442 }
443
444 system_matrix_out.add(index_within_set(
445 index_set_union.nth_index_in_set(row)),
446 temp_indices,
447 temp_values);
448 }
449
450 system_matrix_out.compress(VectorOperation::add);
451 }
452
453
454
455 template <typename SparseMatrixType,
456 typename SparsityPatternType,
457 typename SparseMatrixType2,
458 typename SparsityPatternType2>
459 void
460 restrict_to_serial_sparse_matrix(const SparseMatrixType & system_matrix,
461 const SparsityPatternType &sparsity_pattern,
462 const IndexSet & requested_is,
463 SparseMatrixType2 & system_matrix_out,
464 SparsityPatternType2 &sparsity_pattern_out)
465 {
467 sparsity_pattern,
468 requested_is,
469 IndexSet(), // simply pass empty index set
470 system_matrix_out,
471 sparsity_pattern_out);
472 }
473
474
475
476 template <typename SparseMatrixType,
477 typename SparsityPatternType,
478 typename Number>
479 void
481 const SparseMatrixType & system_matrix,
482 const SparsityPatternType & sparsity_pattern,
483 const std::vector<std::vector<types::global_dof_index>> &indices,
484 std::vector<FullMatrix<Number>> & blocks)
485 {
486 // 0) determine which rows are locally owned and which ones are remote
487 const auto local_size = internal::get_local_size(system_matrix);
488 const auto prefix_sum = internal::compute_prefix_sum(
489 local_size, internal::get_mpi_communicator(system_matrix));
490 IndexSet locally_owned_dofs(std::get<1>(prefix_sum));
491 locally_owned_dofs.add_range(std::get<0>(prefix_sum),
492 std::get<0>(prefix_sum) + local_size);
493
494 std::vector<::types::global_dof_index> ghost_indices_vector;
495
496 for (const auto &i : indices)
497 ghost_indices_vector.insert(ghost_indices_vector.end(),
498 i.begin(),
499 i.end());
500
501 std::sort(ghost_indices_vector.begin(), ghost_indices_vector.end());
502 ghost_indices_vector.erase(std::unique(ghost_indices_vector.begin(),
503 ghost_indices_vector.end()),
504 ghost_indices_vector.end());
505
506
507 IndexSet locally_active_dofs(std::get<1>(prefix_sum));
508 locally_active_dofs.add_indices(ghost_indices_vector.begin(),
509 ghost_indices_vector.end());
510
511 locally_active_dofs.subtract_set(locally_owned_dofs);
512
513 // 1) collect remote rows of sparse matrix
514 const auto locally_relevant_matrix_entries =
515 internal::extract_remote_rows<Number>(system_matrix,
516 sparsity_pattern,
517 locally_active_dofs,
518 internal::get_mpi_communicator(
519 system_matrix));
520
521
522 // 2) loop over all cells and "revert" assembly
523 blocks.clear();
524 blocks.resize(indices.size());
525
526 for (unsigned int c = 0; c < indices.size(); ++c)
527 {
528 if (indices[c].size() == 0)
529 continue;
530
531 const auto &local_dof_indices = indices[c];
532 auto & cell_matrix = blocks[c];
533
534 // allocate memory
535 const unsigned int dofs_per_cell = indices[c].size();
536
537 cell_matrix = FullMatrix<Number>(dofs_per_cell, dofs_per_cell);
538
539 // loop over all entries of the restricted element matrix and
540 // do different things if rows are locally owned or not and
541 // if column entries of that row exist or not
542 for (unsigned int i = 0; i < dofs_per_cell; ++i)
543 for (unsigned int j = 0; j < dofs_per_cell; ++j)
544 {
545 if (locally_owned_dofs.is_element(
546 local_dof_indices[i])) // row is local
547 {
548 cell_matrix(i, j) =
549 sparsity_pattern.exists(local_dof_indices[i],
550 local_dof_indices[j]) ?
551 system_matrix(local_dof_indices[i],
552 local_dof_indices[j]) :
553 0;
554 }
555 else // row is ghost
556 {
557 Assert(locally_active_dofs.is_element(local_dof_indices[i]),
559
560 const auto &row_entries =
561 locally_relevant_matrix_entries[locally_active_dofs
563 local_dof_indices[i])];
564
565 const auto ptr =
566 std::lower_bound(row_entries.begin(),
567 row_entries.end(),
568 std::pair<types::global_dof_index, Number>{
569 local_dof_indices[j], /*dummy*/ 0.0},
570 [](const auto a, const auto b) {
571 return a.first < b.first;
572 });
573
574 if (ptr != row_entries.end() &&
575 local_dof_indices[j] == ptr->first)
576 cell_matrix(i, j) = ptr->second;
577 else
578 cell_matrix(i, j) = 0.0;
579 }
580 }
581 }
582 }
583
584
585
586 template <int dim,
587 int spacedim,
588 typename SparseMatrixType,
589 typename SparsityPatternType,
590 typename Number>
591 void
592 restrict_to_cells(const SparseMatrixType & system_matrix,
593 const SparsityPatternType & sparsity_pattern,
594 const DoFHandler<dim, spacedim> &dof_handler,
595 std::vector<FullMatrix<Number>> &blocks)
596 {
597 std::vector<std::vector<types::global_dof_index>> all_dof_indices;
598 all_dof_indices.resize(dof_handler.get_triangulation().n_active_cells());
599
600 for (const auto &cell : dof_handler.active_cell_iterators())
601 {
602 if (cell->is_locally_owned() == false)
603 continue;
604
605 auto &local_dof_indices = all_dof_indices[cell->active_cell_index()];
606 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
607 cell->get_dof_indices(local_dof_indices);
608 }
609
610 restrict_to_full_matrices(system_matrix,
611 sparsity_pattern,
612 all_dof_indices,
613 blocks);
614 }
615#endif
616
617} // namespace SparseMatrixTools
618
620
621#endif
const Triangulation< dim, spacedim > & get_triangulation() const
size_type size() const
Definition index_set.h:1661
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1883
size_type n_elements() const
Definition index_set.h:1816
bool is_element(const size_type index) const
Definition index_set.h:1776
void subtract_set(const IndexSet &other)
Definition index_set.cc:268
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1727
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
static ::ExceptionBase & ExcInternalError()
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:75
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)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
T sum(const T &t, const MPI_Comm mpi_communicator)
const MPI_Datatype mpi_type_id_for_type
Definition mpi.h:1732
unsigned int global_dof_index
Definition types.h:82
const MPI_Comm comm