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
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
21namespace
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
59namespace 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),
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),
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
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 {
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)
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
BlockType & block(const unsigned int row, const unsigned int column)
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 &)
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 &)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void petsc_increment_state_counter(Vec v)
const MPI_Comm comm