Reference documentation for deal.II version GIT relicensing-439-g5fda5c893d 2024-04-20 06: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\}}\)
Loading...
Searching...
No Matches
petsc_matrix_free.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 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
15
17
18#ifdef DEAL_II_WITH_PETSC
19
22
24
25namespace PETScWrappers
26{
28 {
29 const int m = 0;
30 do_reinit(MPI_COMM_SELF, m, m, m, m);
31 }
32
33
34
35 MatrixFree::MatrixFree(const MPI_Comm communicator,
36 const unsigned int m,
37 const unsigned int n,
38 const unsigned int local_rows,
39 const unsigned int local_columns)
40 {
41 do_reinit(communicator, m, n, local_rows, local_columns);
42 }
43
44
45
47 const MPI_Comm communicator,
48 const unsigned int m,
49 const unsigned int n,
50 const std::vector<unsigned int> &local_rows_per_process,
51 const std::vector<unsigned int> &local_columns_per_process,
52 const unsigned int this_process)
53 {
54 Assert(local_rows_per_process.size() == local_columns_per_process.size(),
55 ExcDimensionMismatch(local_rows_per_process.size(),
56 local_columns_per_process.size()));
57 Assert(this_process < local_rows_per_process.size(), ExcInternalError());
58
59 do_reinit(communicator,
60 m,
61 n,
62 local_rows_per_process[this_process],
63 local_columns_per_process[this_process]);
64 }
65
66
67
68 MatrixFree::MatrixFree(const unsigned int m,
69 const unsigned int n,
70 const unsigned int local_rows,
71 const unsigned int local_columns)
72 {
73 do_reinit(MPI_COMM_WORLD, m, n, local_rows, local_columns);
74 }
75
76
77
79 const unsigned int m,
80 const unsigned int n,
81 const std::vector<unsigned int> &local_rows_per_process,
82 const std::vector<unsigned int> &local_columns_per_process,
83 const unsigned int this_process)
84 {
85 Assert(local_rows_per_process.size() == local_columns_per_process.size(),
86 ExcDimensionMismatch(local_rows_per_process.size(),
87 local_columns_per_process.size()));
88 Assert(this_process < local_rows_per_process.size(), ExcInternalError());
89
90 do_reinit(MPI_COMM_WORLD,
91 m,
92 n,
93 local_rows_per_process[this_process],
94 local_columns_per_process[this_process]);
95 }
96
97
98
99 void
100 MatrixFree::reinit(const MPI_Comm communicator,
101 const unsigned int m,
102 const unsigned int n,
103 const unsigned int local_rows,
104 const unsigned int local_columns)
105 {
106 // destroy the matrix and generate a new one
107 const PetscErrorCode ierr = MatDestroy(&matrix);
108 AssertThrow(ierr == 0, ExcPETScError(ierr));
109
110 do_reinit(communicator, m, n, local_rows, local_columns);
111 }
112
113
114
115 void
116 MatrixFree::reinit(const MPI_Comm communicator,
117 const unsigned int m,
118 const unsigned int n,
119 const std::vector<unsigned int> &local_rows_per_process,
120 const std::vector<unsigned int> &local_columns_per_process,
121 const unsigned int this_process)
122 {
123 Assert(local_rows_per_process.size() == local_columns_per_process.size(),
124 ExcDimensionMismatch(local_rows_per_process.size(),
125 local_columns_per_process.size()));
126 Assert(this_process < local_rows_per_process.size(), ExcInternalError());
127
128 const PetscErrorCode ierr = MatDestroy(&matrix);
129 AssertThrow(ierr != 0, ExcPETScError(ierr));
130
131 do_reinit(communicator,
132 m,
133 n,
134 local_rows_per_process[this_process],
135 local_columns_per_process[this_process]);
136 }
137
138
139
140 void
141 MatrixFree::reinit(const unsigned int m,
142 const unsigned int n,
143 const unsigned int local_rows,
144 const unsigned int local_columns)
145 {
146 reinit(this->get_mpi_communicator(), m, n, local_rows, local_columns);
147 }
148
149
150
151 void
152 MatrixFree::reinit(const unsigned int m,
153 const unsigned int n,
154 const std::vector<unsigned int> &local_rows_per_process,
155 const std::vector<unsigned int> &local_columns_per_process,
156 const unsigned int this_process)
157 {
159 m,
160 n,
161 local_rows_per_process,
162 local_columns_per_process,
163 this_process);
164 }
165
166
167
168 void
170 {
171 const PetscErrorCode ierr = MatDestroy(&matrix);
172 AssertThrow(ierr == 0, ExcPETScError(ierr));
173
174 const int m = 0;
175 do_reinit(MPI_COMM_SELF, m, m, m, m);
176 }
177
178
179
180 void
181 MatrixFree::vmult(Vec &dst, const Vec &src) const
182 {
183 // VectorBase permits us to manipulate, but not own, a Vec
186
187 // This is implemented by derived classes
188 vmult(y, x);
189 }
190
191
192
193 int
194 MatrixFree::matrix_free_mult(Mat A, Vec src, Vec dst)
195 {
196 // create a pointer to this MatrixFree
197 // object and link the given matrix A
198 // to the matrix-vector multiplication
199 // of this MatrixFree object,
200 void *this_object;
201 const PetscErrorCode ierr = MatShellGetContext(A, &this_object);
202 AssertThrow(ierr == 0, ExcPETScError(ierr));
203
204 // call vmult of this object:
205 reinterpret_cast<MatrixFree *>(this_object)->vmult(dst, src);
206
207 return (0);
208 }
209
210
211
212 void
213 MatrixFree::do_reinit(const MPI_Comm communicator,
214 const unsigned int m,
215 const unsigned int n,
216 const unsigned int local_rows,
217 const unsigned int local_columns)
218 {
219 Assert(local_rows <= m, ExcDimensionMismatch(local_rows, m));
220 Assert(local_columns <= n, ExcDimensionMismatch(local_columns, n));
221
222 // create a PETSc MatShell matrix-type
223 // object of dimension m x n and local size
224 // local_rows x local_columns
225 PetscErrorCode ierr = MatCreateShell(communicator,
226 local_rows,
227 local_columns,
228 m,
229 n,
230 static_cast<void *>(this),
231 &matrix);
232 AssertThrow(ierr == 0, ExcPETScError(ierr));
233 // register the MatrixFree::matrix_free_mult function
234 // as the matrix multiplication used by this matrix
235 ierr = MatShellSetOperation(
236 matrix,
237 MATOP_MULT,
238 reinterpret_cast<void (*)()>(
240 AssertThrow(ierr == 0, ExcPETScError(ierr));
241
242 ierr = MatSetFromOptions(matrix);
243 AssertThrow(ierr == 0, ExcPETScError(ierr));
244 }
245} // namespace PETScWrappers
246
247
249
250#endif // DEAL_II_WITH_PETSC
MPI_Comm get_mpi_communicator() const
void reinit(const MPI_Comm communicator, const unsigned int m, const unsigned int n, const unsigned int local_rows, const unsigned int local_columns)
static int matrix_free_mult(Mat A, Vec src, Vec dst)
void do_reinit(const MPI_Comm comm, const unsigned int m, const unsigned int n, const unsigned int local_rows, const unsigned int local_columns)
virtual void vmult(VectorBase &dst, const VectorBase &src) const =0
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define AssertThrow(cond, exc)