deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20:00+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 AssertDimension(row_lengths.size(), m);
155
156 for (const auto &row_length : row_lengths)
157 AssertThrowIntegerConversion(static_cast<PetscInt>(row_length),
158 row_length);
159 const std::vector<PetscInt> int_row_lengths(row_lengths.begin(),
160 row_lengths.end());
161
162 const PetscErrorCode ierr = MatCreateSeqAIJ(
163 PETSC_COMM_SELF, m, n, 0, int_row_lengths.data(), &matrix);
164 AssertThrow(ierr == 0, ExcPETScError(ierr));
165
166 // set symmetric flag, if so requested
167 if (is_symmetric == true)
168 {
169 set_matrix_option(matrix, MAT_SYMMETRIC, PETSC_TRUE);
170 }
171 }
172
173
174
175 template <typename SparsityPatternType>
176 void
177 SparseMatrix::do_reinit(const SparsityPatternType &sparsity_pattern,
178 const bool preset_nonzero_locations)
179 {
180 // If the sparsity pattern's dimensions can be converted to PetscInts then
181 // the rest of the conversions will succeed
182 AssertIntegerConversion(static_cast<PetscInt>(sparsity_pattern.n_rows()),
183 sparsity_pattern.n_rows());
184 AssertIntegerConversion(static_cast<PetscInt>(sparsity_pattern.n_cols()),
185 sparsity_pattern.n_cols());
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]);
217 for (size_type j = 0; j < row_lengths[i]; ++j)
218 {
219 const auto petsc_j =
220 static_cast<PetscInt>(sparsity_pattern.column_number(i, j));
221 row_entries[j] = petsc_j;
222 }
223
224 const auto petsc_i = static_cast<PetscInt>(i);
225 const PetscErrorCode ierr = MatSetValues(matrix,
226 1,
227 &petsc_i,
228 row_lengths[i],
229 row_entries.data(),
230 row_values.data(),
231 INSERT_VALUES);
232 AssertThrow(ierr == 0, ExcPETScError(ierr));
233 }
235
238 }
239 }
240
241 size_t
243 {
244 PetscInt m, n;
245 const PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
246 AssertThrow(ierr == 0, ExcPETScError(ierr));
247
248 return m;
249 }
250
251 size_t
253 {
254 PetscInt m, n;
255 const PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
256 AssertThrow(ierr == 0, ExcPETScError(ierr));
257
258 return n;
259 }
260
261 void
263 const SparseMatrix &B,
264 const MPI::Vector &V) const
265 {
266 // Simply forward to the protected member function of the base class
267 // that takes abstract matrix and vector arguments (to which the compiler
268 // automatically casts the arguments).
269 MatrixBase::mmult(C, B, V);
270 }
271
272 void
274 const SparseMatrix &B,
275 const MPI::Vector &V) const
276 {
277 // Simply forward to the protected member function of the base class
278 // that takes abstract matrix and vector arguments (to which the compiler
279 // automatically casts the arguments).
280 MatrixBase::Tmmult(C, B, V);
281 }
282
283# ifndef DOXYGEN
284 // Explicit instantiations
285 //
286 template SparseMatrix::SparseMatrix(const SparsityPattern &, const bool);
288 const bool);
289
290 template void
291 SparseMatrix::reinit(const SparsityPattern &, const bool);
292 template void
293 SparseMatrix::reinit(const DynamicSparsityPattern &, const bool);
294
295 template void
296 SparseMatrix::do_reinit(const SparsityPattern &, const bool);
297 template void
299# endif
300} // namespace PETScWrappers
301
302
304
305#endif // DEAL_II_WITH_PETSC
size_type row_length(const size_type row) const
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
#define AssertIntegerConversion(index1, index2)
#define AssertThrowIntegerConversion(index1, index2)
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)