Reference documentation for deal.II version 9.5.0
\(\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// 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
25
27
28namespace 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
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
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)
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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)