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