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