Reference documentation for deal.II version GIT c00406db16 2023-09-28 18:00: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_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 
21 # include <deal.II/lac/exceptions.h>
23 
25 
26 namespace 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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
static const char A