Reference documentation for deal.II version GIT relicensing-245-g36f19064f7 2024-03-29 07:20: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\}}\)
Loading...
Searching...
No Matches
petsc_sparse_matrix.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
17#ifdef DEAL_II_WITH_PETSC
18
24
26
27namespace PETScWrappers
28{
30 {
31 const int m = 0, n = 0, n_nonzero_per_row = 0;
32 const PetscErrorCode ierr = MatCreateSeqAIJ(
33 PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix);
34 AssertThrow(ierr == 0, ExcPETScError(ierr));
35 }
36
38 : MatrixBase(A)
39 {}
40
42 const size_type n,
43 const size_type n_nonzero_per_row,
44 const bool is_symmetric)
45 {
46 do_reinit(m, n, n_nonzero_per_row, is_symmetric);
47 }
48
49
50
52 const size_type n,
53 const std::vector<size_type> &row_lengths,
54 const bool is_symmetric)
55 {
56 do_reinit(m, n, row_lengths, is_symmetric);
57 }
58
59
60
61 template <typename SparsityPatternType>
62 SparseMatrix::SparseMatrix(const SparsityPatternType &sparsity_pattern,
63 const bool preset_nonzero_locations)
64 {
65 do_reinit(sparsity_pattern, preset_nonzero_locations);
66 }
67
68
69
72 {
74 return *this;
75 }
76
77
78
79 void
81 const size_type n,
82 const size_type n_nonzero_per_row,
83 const bool is_symmetric)
84 {
85 // get rid of old matrix and generate a
86 // new one
87 const PetscErrorCode ierr = MatDestroy(&matrix);
88 AssertThrow(ierr == 0, ExcPETScError(ierr));
89
90 do_reinit(m, n, n_nonzero_per_row, is_symmetric);
91 }
92
93
94
95 void
97 const size_type n,
98 const std::vector<size_type> &row_lengths,
99 const bool is_symmetric)
100 {
101 // get rid of old matrix and generate a
102 // new one
103 const PetscErrorCode ierr = MatDestroy(&matrix);
104 AssertThrow(ierr == 0, ExcPETScError(ierr));
105
106 do_reinit(m, n, row_lengths, is_symmetric);
107 }
108
109
110
111 template <typename SparsityPatternType>
112 void
113 SparseMatrix::reinit(const SparsityPatternType &sparsity_pattern,
114 const bool preset_nonzero_locations)
115 {
116 // get rid of old matrix and generate a
117 // new one
118 const PetscErrorCode ierr = MatDestroy(&matrix);
119 AssertThrow(ierr == 0, ExcPETScError(ierr));
120
121 do_reinit(sparsity_pattern, preset_nonzero_locations);
122 }
123
124
125
126 void
128 const size_type n,
129 const size_type n_nonzero_per_row,
130 const bool is_symmetric)
131 {
132 // use the call sequence indicating only
133 // a maximal number of elements per row
134 // for all rows globally
135 const PetscErrorCode ierr = MatCreateSeqAIJ(
136 PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix);
137 AssertThrow(ierr == 0, ExcPETScError(ierr));
138
139 // set symmetric flag, if so requested
140 if (is_symmetric == true)
141 {
142 set_matrix_option(matrix, MAT_SYMMETRIC, PETSC_TRUE);
143 }
144 }
145
146
147
148 void
150 const size_type n,
151 const std::vector<size_type> &row_lengths,
152 const bool is_symmetric)
153 {
154 Assert(row_lengths.size() == m,
155 ExcDimensionMismatch(row_lengths.size(), m));
156
157 // use the call sequence indicating a
158 // maximal number of elements for each
159 // row individually. annoyingly, we
160 // always use unsigned ints for cases
161 // like this, while PETSc wants to see
162 // signed integers. so we have to
163 // convert, unless we want to play dirty
164 // tricks with conversions of pointers
165 const std::vector<PetscInt> int_row_lengths(row_lengths.begin(),
166 row_lengths.end());
167
168 const PetscErrorCode ierr = MatCreateSeqAIJ(
169 PETSC_COMM_SELF, m, n, 0, int_row_lengths.data(), &matrix);
170 AssertThrow(ierr == 0, ExcPETScError(ierr));
171
172 // set symmetric flag, if so requested
173 if (is_symmetric == true)
174 {
175 set_matrix_option(matrix, MAT_SYMMETRIC, PETSC_TRUE);
176 }
177 }
178
179
180
181 template <typename SparsityPatternType>
182 void
183 SparseMatrix::do_reinit(const SparsityPatternType &sparsity_pattern,
184 const bool preset_nonzero_locations)
185 {
186 std::vector<size_type> row_lengths(sparsity_pattern.n_rows());
187 for (size_type i = 0; i < sparsity_pattern.n_rows(); ++i)
188 row_lengths[i] = sparsity_pattern.row_length(i);
189
190 do_reinit(sparsity_pattern.n_rows(),
191 sparsity_pattern.n_cols(),
192 row_lengths,
193 false);
194
195 // next preset the exact given matrix
196 // entries with zeros, if the user
197 // requested so. this doesn't avoid any
198 // memory allocations, but it at least
199 // avoids some searches later on. the
200 // key here is that we can use the
201 // matrix set routines that set an
202 // entire row at once, not a single
203 // entry at a time
204 //
205 // for the usefulness of this option
206 // read the documentation of this
207 // class.
208 if (preset_nonzero_locations == true)
209 {
210 std::vector<PetscInt> row_entries;
211 std::vector<PetscScalar> row_values;
212 for (size_type i = 0; i < sparsity_pattern.n_rows(); ++i)
213 {
214 row_entries.resize(row_lengths[i]);
215 row_values.resize(row_lengths[i], 0.0);
216 for (size_type j = 0; j < row_lengths[i]; ++j)
217 row_entries[j] = sparsity_pattern.column_number(i, j);
218
219 const PetscInt int_row = i;
220 const PetscErrorCode ierr = MatSetValues(matrix,
221 1,
222 &int_row,
223 row_lengths[i],
224 row_entries.data(),
225 row_values.data(),
226 INSERT_VALUES);
227 AssertThrow(ierr == 0, ExcPETScError(ierr));
228 }
230
233 }
234 }
235
236 size_t
238 {
239 PetscInt m, n;
240 const PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
241 AssertThrow(ierr == 0, ExcPETScError(ierr));
242
243 return m;
244 }
245
246 size_t
248 {
249 PetscInt m, n;
250 const PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
251 AssertThrow(ierr == 0, ExcPETScError(ierr));
252
253 return n;
254 }
255
256 void
258 const SparseMatrix &B,
259 const MPI::Vector &V) const
260 {
261 // Simply forward to the protected member function of the base class
262 // that takes abstract matrix and vector arguments (to which the compiler
263 // automatically casts the arguments).
264 MatrixBase::mmult(C, B, V);
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::Tmmult(C, B, V);
276 }
277
278# ifndef DOXYGEN
279 // Explicit instantiations
280 //
281 template SparseMatrix::SparseMatrix(const SparsityPattern &, const bool);
283 const bool);
284
285 template void
286 SparseMatrix::reinit(const SparsityPattern &, const bool);
287 template void
288 SparseMatrix::reinit(const DynamicSparsityPattern &, const bool);
289
290 template void
291 SparseMatrix::do_reinit(const SparsityPattern &, const bool);
292 template void
294# endif
295} // namespace PETScWrappers
296
297
299
300#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)
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define AssertThrow(cond, exc)
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)