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