Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
petsc_sparse_matrix.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2018 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 #include <deal.II/lac/petsc_sparse_matrix.h>
17 
18 #ifdef DEAL_II_WITH_PETSC
19 
20 # include <deal.II/lac/dynamic_sparsity_pattern.h>
21 # include <deal.II/lac/exceptions.h>
22 # include <deal.II/lac/petsc_compatibility.h>
23 # include <deal.II/lac/petsc_vector_base.h>
24 # include <deal.II/lac/sparsity_pattern.h>
25 
26 DEAL_II_NAMESPACE_OPEN
27 
28 namespace PETScWrappers
29 {
31  {
32  const int m = 0, n = 0, n_nonzero_per_row = 0;
33  const PetscErrorCode ierr = MatCreateSeqAIJ(
34  PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix);
35  AssertThrow(ierr == 0, ExcPETScError(ierr));
36  }
37 
38 
39 
41  const size_type n,
42  const size_type n_nonzero_per_row,
43  const bool is_symmetric)
44  {
45  do_reinit(m, n, n_nonzero_per_row, is_symmetric);
46  }
47 
48 
49 
51  const size_type n,
52  const std::vector<size_type> &row_lengths,
53  const bool is_symmetric)
54  {
55  do_reinit(m, n, row_lengths, is_symmetric);
56  }
57 
58 
59 
60  template <typename SparsityPatternType>
61  SparseMatrix::SparseMatrix(const SparsityPatternType &sparsity_pattern,
62  const bool preset_nonzero_locations)
63  {
64  do_reinit(sparsity_pattern, preset_nonzero_locations);
65  }
66 
67 
68 
69  SparseMatrix &
70  SparseMatrix::operator=(const double d)
71  {
73  return *this;
74  }
75 
76 
77 
78  void
80  const size_type n,
81  const size_type n_nonzero_per_row,
82  const bool is_symmetric)
83  {
84  // get rid of old matrix and generate a
85  // new one
86  const PetscErrorCode ierr = destroy_matrix(matrix);
87  AssertThrow(ierr == 0, ExcPETScError(ierr));
88 
89  do_reinit(m, n, n_nonzero_per_row, is_symmetric);
90  }
91 
92 
93 
94  void
96  const size_type n,
97  const std::vector<size_type> &row_lengths,
98  const bool is_symmetric)
99  {
100  // get rid of old matrix and generate a
101  // new one
102  const PetscErrorCode ierr = destroy_matrix(matrix);
103  AssertThrow(ierr == 0, ExcPETScError(ierr));
104 
105  do_reinit(m, n, row_lengths, is_symmetric);
106  }
107 
108 
109 
110  template <typename SparsityPatternType>
111  void
112  SparseMatrix::reinit(const SparsityPatternType &sparsity_pattern,
113  const bool preset_nonzero_locations)
114  {
115  // get rid of old matrix and generate a
116  // new one
117  const PetscErrorCode ierr = destroy_matrix(matrix);
118  AssertThrow(ierr == 0, ExcPETScError(ierr));
119 
120  do_reinit(sparsity_pattern, preset_nonzero_locations);
121  }
122 
123 
124 
125  const MPI_Comm &
127  {
128  static MPI_Comm comm;
129  const PetscErrorCode ierr =
130  PetscObjectGetComm(reinterpret_cast<PetscObject>(matrix), &comm);
131  AssertThrow(ierr == 0, ExcPETScError(ierr));
132  return comm;
133  }
134 
135 
136 
137  void
139  const size_type n,
140  const size_type n_nonzero_per_row,
141  const bool is_symmetric)
142  {
143  // use the call sequence indicating only
144  // a maximal number of elements per row
145  // for all rows globally
146  const PetscErrorCode ierr = MatCreateSeqAIJ(
147  PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix);
148  AssertThrow(ierr == 0, ExcPETScError(ierr));
149 
150  // set symmetric flag, if so requested
151  if (is_symmetric == true)
152  {
153  set_matrix_option(matrix, MAT_SYMMETRIC, PETSC_TRUE);
154  }
155  }
156 
157 
158 
159  void
161  const size_type n,
162  const std::vector<size_type> &row_lengths,
163  const bool is_symmetric)
164  {
165  Assert(row_lengths.size() == m,
166  ExcDimensionMismatch(row_lengths.size(), m));
167 
168  // use the call sequence indicating a
169  // maximal number of elements for each
170  // row individually. annoyingly, we
171  // always use unsigned ints for cases
172  // like this, while PETSc wants to see
173  // signed integers. so we have to
174  // convert, unless we want to play dirty
175  // tricks with conversions of pointers
176  const std::vector<PetscInt> int_row_lengths(row_lengths.begin(),
177  row_lengths.end());
178 
179  const PetscErrorCode ierr = MatCreateSeqAIJ(
180  PETSC_COMM_SELF, m, n, 0, int_row_lengths.data(), &matrix);
181  AssertThrow(ierr == 0, ExcPETScError(ierr));
182 
183  // set symmetric flag, if so requested
184  if (is_symmetric == true)
185  {
186  set_matrix_option(matrix, MAT_SYMMETRIC, PETSC_TRUE);
187  }
188  }
189 
190 
191 
192  template <typename SparsityPatternType>
193  void
194  SparseMatrix::do_reinit(const SparsityPatternType &sparsity_pattern,
195  const bool preset_nonzero_locations)
196  {
197  std::vector<size_type> row_lengths(sparsity_pattern.n_rows());
198  for (size_type i = 0; i < sparsity_pattern.n_rows(); ++i)
199  row_lengths[i] = sparsity_pattern.row_length(i);
200 
201  do_reinit(sparsity_pattern.n_rows(),
202  sparsity_pattern.n_cols(),
203  row_lengths,
204  false);
205 
206  // next preset the exact given matrix
207  // entries with zeros, if the user
208  // requested so. this doesn't avoid any
209  // memory allocations, but it at least
210  // avoids some searches later on. the
211  // key here is that we can use the
212  // matrix set routines that set an
213  // entire row at once, not a single
214  // entry at a time
215  //
216  // for the usefulness of this option
217  // read the documentation of this
218  // class.
219  if (preset_nonzero_locations == true)
220  {
221  std::vector<PetscInt> row_entries;
222  std::vector<PetscScalar> row_values;
223  for (size_type i = 0; i < sparsity_pattern.n_rows(); ++i)
224  {
225  row_entries.resize(row_lengths[i]);
226  row_values.resize(row_lengths[i], 0.0);
227  for (size_type j = 0; j < row_lengths[i]; ++j)
228  row_entries[j] = sparsity_pattern.column_number(i, j);
229 
230  const PetscInt int_row = i;
231  const PetscErrorCode ierr = MatSetValues(matrix,
232  1,
233  &int_row,
234  row_lengths[i],
235  row_entries.data(),
236  row_values.data(),
237  INSERT_VALUES);
238  AssertThrow(ierr == 0, ExcPETScError(ierr));
239  }
241 
244  }
245  }
246 
247  size_t
249  {
250  PetscInt m, n;
251  const PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
252  AssertThrow(ierr == 0, ExcPETScError(ierr));
253 
254  return m;
255  }
256 
257  size_t
259  {
260  PetscInt m, n;
261  const PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
262  AssertThrow(ierr == 0, ExcPETScError(ierr));
263 
264  return n;
265  }
266 
267  void
269  const SparseMatrix &B,
270  const MPI::Vector & V) const
271  {
272  // Simply forward to the protected member function of the base class
273  // that takes abstract matrix and vector arguments (to which the compiler
274  // automatically casts the arguments).
275  MatrixBase::mmult(C, B, V);
276  }
277 
278  void
280  const SparseMatrix &B,
281  const MPI::Vector & V) const
282  {
283  // Simply forward to the protected member function of the base class
284  // that takes abstract matrix and vector arguments (to which the compiler
285  // automatically casts the arguments).
286  MatrixBase::Tmmult(C, B, V);
287  }
288 
289  // Explicit instantiations
290  //
291  template SparseMatrix::SparseMatrix(const SparsityPattern &, const bool);
293  const bool);
294 
295  template void
296  SparseMatrix::reinit(const SparsityPattern &, const bool);
297  template void
298  SparseMatrix::reinit(const DynamicSparsityPattern &, const bool);
299 
300  template void
301  SparseMatrix::do_reinit(const SparsityPattern &, const bool);
302  template void
303  SparseMatrix::do_reinit(const DynamicSparsityPattern &, const bool);
304 } // namespace PETScWrappers
305 
306 
307 DEAL_II_NAMESPACE_CLOSE
308 
309 #endif // DEAL_II_WITH_PETSC
PetscErrorCode destroy_matrix(Mat &matrix)
void Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
virtual const MPI_Comm & get_mpi_communicator() const override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
void set_matrix_option(Mat &matrix, const MatOption option_name, const PetscBool option_value=PETSC_FALSE)
PetscBool is_symmetric(const double tolerance=1.e-12)
void compress(const VectorOperation::values operation)
MatrixBase & operator=(const MatrixBase &)=delete
SparseMatrix & operator=(const double d)
void reinit(const size_type m, const size_type n, const size_type n_nonzero_per_row, const bool is_symmetric=false)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
void set_keep_zero_rows(Mat &matrix)
types::global_dof_index size_type
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
void do_reinit(const size_type m, const size_type n, const size_type n_nonzero_per_row, const bool is_symmetric=false)
void close_matrix(Mat &matrix)