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