Reference documentation for deal.II version GIT f6a5d312c9 2023-10-04 08:50: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\}}\)
petsc_parallel_block_sparse_matrix.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 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 
18 
19 #ifdef DEAL_II_WITH_PETSC
20 
21 namespace
22 {
23  // A dummy utility routine to create an empty matrix in case we import
24  // a MATNEST with NULL blocks
25  static Mat
26  create_dummy_mat(MPI_Comm comm,
27  PetscInt lr,
28  PetscInt gr,
29  PetscInt lc,
30  PetscInt gc)
31  {
32  Mat dummy;
33  PetscErrorCode ierr;
34 
35  ierr = MatCreate(comm, &dummy);
36  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
37  ierr = MatSetSizes(dummy, lr, lc, gr, gc);
38  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
39  ierr = MatSetType(dummy, MATAIJ);
40  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
41  ierr = MatSeqAIJSetPreallocation(dummy, 0, nullptr);
42  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
43  ierr = MatMPIAIJSetPreallocation(dummy, 0, nullptr, 0, nullptr);
44  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
45  ierr = MatSetUp(dummy);
46  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
47  ierr = MatSetOption(dummy, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
48  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
49  ierr = MatAssemblyBegin(dummy, MAT_FINAL_ASSEMBLY);
50  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
51  ierr = MatAssemblyEnd(dummy, MAT_FINAL_ASSEMBLY);
52  AssertThrow(ierr == 0, ::ExcPETScError(ierr));
53  return dummy;
54  }
55 } // namespace
56 
58 
59 namespace PETScWrappers
60 {
61  namespace MPI
62  {
65  {
67 
68  return *this;
69  }
70 
71 
72 
74  {
75  PetscErrorCode ierr = MatDestroy(&petsc_nest_matrix);
76  (void)ierr;
77  AssertNothrow(ierr == 0, ExcPETScError(ierr));
78  }
79 
80 
81 
82 # ifndef DOXYGEN
83  void
84  BlockSparseMatrix::reinit(const size_type n_block_rows,
85  const size_type n_block_columns)
86  {
87  // first delete previous content of
88  // the subobjects array
89  clear();
90 
91  // then resize. set sizes of blocks to
92  // zero. user will later have to call
93  // collect_sizes for this
94  this->sub_objects.reinit(n_block_rows, n_block_columns);
95  this->row_block_indices.reinit(n_block_rows, 0);
96  this->column_block_indices.reinit(n_block_columns, 0);
97 
98  // and reinitialize the blocks
99  for (size_type r = 0; r < this->n_block_rows(); ++r)
100  for (size_type c = 0; c < this->n_block_cols(); ++c)
101  {
102  BlockType *p = new BlockType();
103  this->sub_objects[r][c] = p;
104  }
105  }
106 # endif
107 
108 
109 
110  void
111  BlockSparseMatrix::reinit(const std::vector<IndexSet> &rows,
112  const std::vector<IndexSet> &cols,
113  const BlockDynamicSparsityPattern &bdsp,
114  const MPI_Comm com)
115  {
116  Assert(rows.size() == bdsp.n_block_rows(), ExcMessage("invalid size"));
117  Assert(cols.size() == bdsp.n_block_cols(), ExcMessage("invalid size"));
118 
119 
120  clear();
121  this->sub_objects.reinit(bdsp.n_block_rows(), bdsp.n_block_cols());
122 
123  std::vector<types::global_dof_index> row_sizes;
124  for (unsigned int r = 0; r < bdsp.n_block_rows(); ++r)
125  row_sizes.push_back(bdsp.block(r, 0).n_rows());
126  this->row_block_indices.reinit(row_sizes);
127 
128  std::vector<types::global_dof_index> col_sizes;
129  for (unsigned int c = 0; c < bdsp.n_block_cols(); ++c)
130  col_sizes.push_back(bdsp.block(0, c).n_cols());
131  this->column_block_indices.reinit(col_sizes);
132 
133  for (unsigned int r = 0; r < this->n_block_rows(); ++r)
134  for (unsigned int c = 0; c < this->n_block_cols(); ++c)
135  {
136  Assert(rows[r].size() == bdsp.block(r, c).n_rows(),
137  ExcMessage("invalid size"));
138  Assert(cols[c].size() == bdsp.block(r, c).n_cols(),
139  ExcMessage("invalid size"));
140 
141  BlockType *p = new BlockType();
142  p->reinit(rows[r], cols[c], bdsp.block(r, c), com);
143  this->sub_objects[r][c] = p;
144  }
145 
146  this->collect_sizes();
147  }
148 
149  void
150  BlockSparseMatrix::reinit(const std::vector<IndexSet> &sizes,
151  const BlockDynamicSparsityPattern &bdsp,
152  const MPI_Comm com)
153  {
154  reinit(sizes, sizes, bdsp, com);
155  }
156 
157 
158 
159  void
161  {
162  auto m = this->n_block_rows();
163  auto n = this->n_block_cols();
164  PetscErrorCode ierr;
165 
166  // Create empty matrices if needed
167  // This is needed by the base class
168  // not by MATNEST
169  std::vector<size_type> row_sizes(m, size_type(-1));
170  std::vector<size_type> col_sizes(n, size_type(-1));
171  std::vector<size_type> row_local_sizes(m, size_type(-1));
172  std::vector<size_type> col_local_sizes(n, size_type(-1));
173  MPI_Comm comm = MPI_COMM_NULL;
174  for (size_type r = 0; r < m; r++)
175  {
176  for (size_type c = 0; c < n; c++)
177  {
178  if (this->sub_objects[r][c])
179  {
180  comm = this->sub_objects[r][c]->get_mpi_communicator();
181  row_sizes[r] = this->sub_objects[r][c]->m();
182  col_sizes[c] = this->sub_objects[r][c]->n();
183  row_local_sizes[r] = this->sub_objects[r][c]->local_size();
184  col_local_sizes[c] =
185  this->sub_objects[r][c]->local_domain_size();
186  }
187  }
188  }
189  for (size_type r = 0; r < m; r++)
190  {
191  for (size_type c = 0; c < n; c++)
192  {
193  if (!this->sub_objects[r][c])
194  {
195  Assert(
196  row_sizes[r] != size_type(-1),
197  ExcMessage(
198  "When passing empty sub-blocks of a block matrix, you need to make "
199  "sure that at least one block in each block row and block column is "
200  "non-empty. However, block row " +
201  std::to_string(r) +
202  " is completely empty "
203  "and so it is not possible to determine how many rows it should have."));
204  Assert(
205  col_sizes[c] != size_type(-1),
206  ExcMessage(
207  "When passing empty sub-blocks of a block matrix, you need to make "
208  "sure that at least one block in each block row and block column is "
209  "non-empty. However, block column " +
210  std::to_string(c) +
211  " is completely empty "
212  "and so it is not possible to determine how many columns it should have."));
213  Mat dummy = ::create_dummy_mat(
214  comm,
215  static_cast<PetscInt>(row_local_sizes[r]),
216  static_cast<PetscInt>(row_sizes[r]),
217  static_cast<PetscInt>(col_local_sizes[c]),
218  static_cast<PetscInt>(col_sizes[c]));
219  this->sub_objects[r][c] = new BlockType(dummy);
220 
221  // the new object got a reference on dummy, we can safely
222  // call destroy here
223  ierr = MatDestroy(&dummy);
224  AssertThrow(ierr == 0, ExcPETScError(ierr));
225  }
226  }
227  }
228  }
229 
230 
231  void
233  {
236  this->setup_nest_mat();
237  }
238 
239  void
241  {
242  auto m = this->n_block_rows();
243  auto n = this->n_block_cols();
244  PetscErrorCode ierr;
245 
246  MPI_Comm comm = PETSC_COMM_SELF;
247 
248  ierr = MatDestroy(&petsc_nest_matrix);
249  AssertThrow(ierr == 0, ExcPETScError(ierr));
250  std::vector<Mat> psub_objects(m * n);
251  for (unsigned int r = 0; r < m; r++)
252  for (unsigned int c = 0; c < n; c++)
253  {
254  comm = this->sub_objects[r][c]->get_mpi_communicator();
255  psub_objects[r * n + c] = this->sub_objects[r][c]->petsc_matrix();
256  }
257  ierr = MatCreateNest(
258  comm, m, nullptr, n, nullptr, psub_objects.data(), &petsc_nest_matrix);
259  AssertThrow(ierr == 0, ExcPETScError(ierr));
260 
261  ierr = MatNestSetVecType(petsc_nest_matrix, VECNEST);
262  AssertThrow(ierr == 0, ExcPETScError(ierr));
263  }
264 
265 
266 
267  void
269  {
270  BaseClass::compress(operation);
272  }
273 
274 
275 
276  std::vector<IndexSet>
278  {
279  std::vector<IndexSet> index_sets;
280 
281  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
282  index_sets.push_back(this->block(0, i).locally_owned_domain_indices());
283 
284  return index_sets;
285  }
286 
287 
288 
289  std::vector<IndexSet>
291  {
292  std::vector<IndexSet> index_sets;
293 
294  for (unsigned int i = 0; i < this->n_block_rows(); ++i)
295  index_sets.push_back(this->block(i, 0).locally_owned_range_indices());
296 
297  return index_sets;
298  }
299 
300 
301 
302  std::uint64_t
304  {
305  std::uint64_t n_nonzero = 0;
306  for (size_type rows = 0; rows < this->n_block_rows(); ++rows)
307  for (size_type cols = 0; cols < this->n_block_cols(); ++cols)
308  n_nonzero += this->block(rows, cols).n_nonzero_elements();
309 
310  return n_nonzero;
311  }
312 
313 
314 
315  MPI_Comm
317  {
318  return PetscObjectComm(reinterpret_cast<PetscObject>(petsc_nest_matrix));
319  }
320 
321  BlockSparseMatrix::operator const Mat &() const
322  {
323  return petsc_nest_matrix;
324  }
325 
326 
327 
328  Mat &
330  {
331  return petsc_nest_matrix;
332  }
333 
334  void
336  {
337  clear();
338 
339  PetscBool isnest;
340  PetscInt nr = 1, nc = 1;
341 
342  PetscErrorCode ierr =
343  PetscObjectTypeCompare(reinterpret_cast<PetscObject>(A),
344  MATNEST,
345  &isnest);
346  AssertThrow(ierr == 0, ExcPETScError(ierr));
347  std::vector<Mat> mats;
348  bool need_empty_matrices = false;
349  if (isnest)
350  {
351  ierr = MatNestGetSize(A, &nr, &nc);
352  AssertThrow(ierr == 0, ExcPETScError(ierr));
353  for (PetscInt i = 0; i < nr; ++i)
354  {
355  for (PetscInt j = 0; j < nc; ++j)
356  {
357  Mat sA;
358  ierr = MatNestGetSubMat(A, i, j, &sA);
359  mats.push_back(sA);
360  if (!sA)
361  need_empty_matrices = true;
362  }
363  }
364  }
365  else
366  {
367  mats.push_back(A);
368  }
369 
370  std::vector<size_type> r_block_sizes(nr, 0);
371  std::vector<size_type> c_block_sizes(nc, 0);
372  this->row_block_indices.reinit(r_block_sizes);
373  this->column_block_indices.reinit(c_block_sizes);
374  this->sub_objects.reinit(nr, nc);
375  for (PetscInt i = 0; i < nr; ++i)
376  {
377  for (PetscInt j = 0; j < nc; ++j)
378  {
379  if (mats[i * nc + j])
380  this->sub_objects[i][j] = new BlockType(mats[i * nc + j]);
381  else
382  this->sub_objects[i][j] = nullptr;
383  }
384  }
385  if (need_empty_matrices)
387 
389  if (need_empty_matrices || !isnest)
390  {
391  setup_nest_mat();
392  }
393  else
394  {
395  ierr = PetscObjectReference(reinterpret_cast<PetscObject>(A));
396  AssertThrow(ierr == 0, ExcPETScError(ierr));
397  PetscErrorCode ierr = MatDestroy(&petsc_nest_matrix);
398  AssertThrow(ierr == 0, ExcPETScError(ierr));
400  }
401  }
402 
403  } // namespace MPI
404 } // namespace PETScWrappers
405 
406 
407 
409 
410 #endif
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
BlockType & block(const unsigned int row, const unsigned int column)
unsigned int n_block_rows() const
void compress(VectorOperation::values operation)
Table< 2, SmartPointer< BlockType, BlockMatrixBase< SparseMatrix > > > sub_objects
unsigned int n_block_cols() const
SparsityPatternType & block(const size_type row, const size_type column)
void reinit(const size_type n_block_rows, const size_type n_block_columns)
BlockSparseMatrix & operator=(const BlockSparseMatrix &)
std::vector< IndexSet > locally_owned_domain_indices() const
void compress(VectorOperation::values operation)
std::size_t n_nonzero_elements() const
virtual void reinit(const SparsityPattern &sparsity)
size_type n_rows() const
size_type n_cols() const
Subscriptor & operator=(const Subscriptor &)
Definition: subscriptor.h:291
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1658
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
static const char A
types::global_dof_index size_type
void petsc_increment_state_counter(Vec v)
std::string to_string(const T &t)
Definition: patterns.h:2391
const MPI_Comm comm