AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
MatrixBackend.hpp
1#pragma once
2
3#include <algorithm>
4#include <memory>
5#include <vector>
6
7#include <petscmat.h>
8
9#include <dune/common/timer.hh>
10
11#include <amdis/Output.hpp>
12#include <amdis/common/FakeContainer.hpp>
13#include <amdis/linearalgebra/petsc/MatrixNnzStructure.hpp>
14#include <amdis/linearalgebra/petsc/VectorBackend.hpp>
15#include <amdis/linearalgebra/SymmetryStructure.hpp>
16
17namespace AMDiS
18{
19 template <class> class Constraints;
20
22 template <class DofMap>
24 {
25 template <class> friend class Constraints;
26
27 public:
29 using BaseMatrix = ::Mat;
30
32 using value_type = PetscScalar;
33
35 using size_type = PetscInt;
36
37 public:
39 PetscMatrix(DofMap const& rowDofMap, DofMap const& colDofMap)
40 : rowDofMap_(rowDofMap)
41 , colDofMap_(colDofMap)
42 {}
43
44 // disable copy and move operations
45 PetscMatrix(PetscMatrix const&) = delete;
46 PetscMatrix(PetscMatrix&&) = delete;
47 PetscMatrix& operator=(PetscMatrix const&) = delete;
48 PetscMatrix& operator=(PetscMatrix&&) = delete;
49
51 {
52 destroy();
53 }
54
57 {
58 return matrix_;
59 }
60
62 BaseMatrix const& matrix() const
63 {
64 return matrix_;
65 }
66
68 template <class RowIndex, class ColIndex>
69 void insert(RowIndex const& r, ColIndex const& c, PetscScalar value)
70 {
71 PetscInt row = rowDofMap_.global(r);
72 PetscInt col = colDofMap_.global(c);
73 MatSetValue(matrix_, row, col, value, ADD_VALUES);
74 }
75
77 template <class LocalInd, class LocalMatrix>
78 void scatter(LocalInd const& localInd, LocalMatrix const& localMat)
79 {
80 if (&rowDofMap_ == &colDofMap_) {
81 thread_local std::vector<PetscInt> idx;
82 idx.resize(localInd.size());
83
84 // create a vector of global indices from the local indices using the local-to-global map
85 std::transform(localInd.begin(), localInd.end(), idx.begin(),
86 [this](auto const& mi) { return rowDofMap_.global(mi); });
87
88 MatSetValues(matrix_, idx.size(), idx.data(), idx.size(), idx.data(), localMat.data(), ADD_VALUES);
89 } else {
90 scatter(localInd, localInd, localMat);
91 }
92 }
93
95 template <class RowLocalInd, class ColLocalInd, class LocalMatrix>
96 void scatter(RowLocalInd const& rowLocalInd, ColLocalInd const& colLocalInd, LocalMatrix const& localMat)
97 {
98 thread_local std::vector<PetscInt> ri;
99 thread_local std::vector<PetscInt> ci;
100 ri.resize(rowLocalInd.size());
101 ci.resize(colLocalInd.size());
102
103 // create vectors of global indices from the local indices using the local-to-global map
104 std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
105 [this](auto const& mi) { return rowDofMap_.global(mi); });
106 std::transform(colLocalInd.begin(), colLocalInd.end(), ci.begin(),
107 [this](auto const& mi) { return colDofMap_.global(mi); });
108
109 MatSetValues(matrix_, ri.size(), ri.data(), ci.size(), ci.data(), localMat.data(), ADD_VALUES);
110 }
111
113 template <class RowInd>
114 void zeroRows(std::vector<RowInd> const& rowLocalInd, bool diag)
115 {
116 thread_local std::vector<PetscInt> ri;
117 ri.resize(rowLocalInd.size());
118 std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
119 [this](auto const& mi) { return rowDofMap_.global(mi); });
120
121 MatZeroRows(matrix_, ri.size(), ri.data(), diag ? 1.0 : 0.0, PETSC_NULL, PETSC_NULL);
122 }
123
125 template <class RowInd, class ColInd>
126 void zeroRows(std::vector<RowInd> const& rowLocalInd, std::vector<ColInd> const& colLocalInd, bool diag)
127 {
128 thread_local std::vector<PetscInt> ri;
129 ri.resize(rowLocalInd.size());
130 std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
131 [this](auto const& mi) { return rowDofMap_.global(mi); });
132 MatZeroRows(matrix_, ri.size(), ri.data(), 0.0, PETSC_NULL, PETSC_NULL);
133
134 if (diag) {
135 thread_local std::vector<PetscInt> ci;
136 ci.resize(colLocalInd.size());
137 std::transform(colLocalInd.begin(), colLocalInd.end(), ci.begin(),
138 [this](auto const& mi) { return colDofMap_.global(mi); });
139 PetscBool mat_new_nonzero_allocation_err;
140 MatGetOption(matrix_, MAT_NEW_NONZERO_ALLOCATION_ERR, &mat_new_nonzero_allocation_err);
141 MatSetOption(matrix_, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
142 for (std::size_t i = 0; i < ri.size(); ++i)
143 MatSetValue(matrix_, ri[i], ci[i], 1.0, INSERT_VALUES);
144 MatAssemblyBegin(matrix_, MAT_FINAL_ASSEMBLY);
145 MatAssemblyEnd(matrix_, MAT_FINAL_ASSEMBLY);
146 MatSetOption(matrix_, MAT_NEW_NONZERO_ALLOCATION_ERR, mat_new_nonzero_allocation_err);
147 }
148 }
149
151 template <class RowInd>
152 void zeroRowsColumns(std::vector<RowInd> const& rowLocalInd, bool diag)
153 {
154 thread_local std::vector<PetscInt> ri;
155 ri.resize(rowLocalInd.size());
156 std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
157 [this](auto const& mi) { return rowDofMap_.global(mi); });
158 MatZeroRowsColumns(matrix_, ri.size(), ri.data(), diag ? 1.0 : 0.0, PETSC_NULL, PETSC_NULL);
159 }
160
162 template <class RowInd>
163 void zeroRowsColumns(std::vector<RowInd> const& rowLocalInd, bool diag, PetscVector<DofMap> const& x, PetscVector<DofMap>& b)
164 {
165 thread_local std::vector<PetscInt> ri;
166 ri.resize(rowLocalInd.size());
167 std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
168 [this](auto const& mi) { return rowDofMap_.global(mi); });
169 MatZeroRowsColumns(matrix_, ri.size(), ri.data(), diag ? 1.0 : 0.0, x.vector(), b.vector());
170 }
171
173 template <class RowInd, class ColInd>
174 void zeroRowsColumns(std::vector<RowInd> const& rowLocalInd, std::vector<ColInd> const& colLocalInd, bool diag)
175 {
176 error_exit("Different row and column indices not implemented for PetscMatrix::zeroRowsColumns. Either pass only one index vector, if the row and column indices are the same, or use PetscMatrix::zeroRows that does not eliminate the columns.");
177 }
178
180 template <class RowInd, class ColInd>
181 void zeroRowsColumns(std::vector<RowInd> const& rowLocalInd, std::vector<ColInd> const& colLocalInd, bool diag, PetscVector<DofMap> const& x, PetscVector<DofMap>& b)
182 {
183 error_exit("Different row and column indices not implemented for PetscMatrix::zeroRowsColumns. Either pass only one index vector, if the row and column indices are the same, or use PetscMatrix::zeroRows that does not eliminate the columns.");
184 }
185
187 void init(MatrixNnzStructure const& pattern)
188 {
189 Dune::Timer t;
190
191 // destroy an old matrix if created before
192 destroy();
193 info(3, " destroy old matrix needed {} seconds", t.elapsed());
194 t.reset();
195
196 // create a MATAIJ or MATSEQAIJ matrix with provided sparsity pattern
197 MatCreateAIJ(comm(),
198 rowDofMap_.localSize(), colDofMap_.localSize(),
199 rowDofMap_.globalSize(), colDofMap_.globalSize(),
200 0, pattern.d_nnz().data(), 0, pattern.o_nnz().data(), &matrix_);
201
202 // keep sparsity pattern even if we delete a row / column with e.g. MatZeroRows
203 MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
204
205 // set symmetry properties of the matrix
206 switch (pattern.symmetry()) {
207 case SymmetryStructure::spd:
208 MatSetOption(matrix_, MAT_SPD, PETSC_TRUE);
209 break;
210 case SymmetryStructure::symmetric:
211 MatSetOption(matrix_, MAT_SYMMETRIC, PETSC_TRUE);
212 break;
213 case SymmetryStructure::hermitian:
214 MatSetOption(matrix_, MAT_HERMITIAN, PETSC_TRUE);
215 break;
216 case SymmetryStructure::structurally_symmetric:
217 MatSetOption(matrix_, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
218 break;
219 default:
220 /* do nothing */
221 break;
222 }
223
224 info(3, " create new matrix needed {} seconds", t.elapsed());
225 t.reset();
226
227 initialized_ = true;
228 }
229
231 void init()
232 {
233 MatZeroEntries(matrix_);
234 initialized_ = true;
235 }
236
238 void finish()
239 {
240 Dune::Timer t;
241 MatAssemblyBegin(matrix_, MAT_FINAL_ASSEMBLY);
242 MatAssemblyEnd(matrix_, MAT_FINAL_ASSEMBLY);
243 info(3, " finish matrix assembling needed {} seconds", t.elapsed());
244 }
245
247 std::size_t nnz() const
248 {
249 MatInfo info;
250 MatGetInfo(matrix_, MAT_LOCAL, &info);
251 return std::size_t(info.nz_used);
252 }
253
254 // An MPI Communicator or PETSC_COMM_SELF
255 MPI_Comm comm() const
256 {
257 return rowDofMap_.comm();
258 }
259
260 private:
261 // Destroy the matrix if initialized before.
262 void destroy()
263 {
264 if (initialized_)
265 MatDestroy(&matrix_);
266 }
267
268 private:
269 // The local-to-global maps
270 DofMap const& rowDofMap_;
271 DofMap const& colDofMap_;
272
274 Mat matrix_;
275
276 bool initialized_ = false;
277 };
278
279} // end namespace AMDiS
Sparsity pattern used to create PETSc matrices.
Definition: MatrixNnzStructure.hpp:18
SymmetryStructure symmetry() const
Symmetry of the matrix entries.
Definition: MatrixNnzStructure.hpp:41
std::vector< PetscInt > const & o_nnz() const
Return Number of nonzeros in the off-diagonal part (overlap part)
Definition: MatrixNnzStructure.hpp:35
The basic container that stores a base matrix.
Definition: MatrixBackend.hpp:24
void init()
Reuse the matrix pattern and set all entries to zero.
Definition: MatrixBackend.hpp:231
std::size_t nnz() const
Return the local number of nonzeros in the matrix.
Definition: MatrixBackend.hpp:247
void insert(RowIndex const &r, ColIndex const &c, PetscScalar value)
Insert a single value into the matrix.
Definition: MatrixBackend.hpp:69
BaseMatrix const & matrix() const
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:62
void scatter(LocalInd const &localInd, LocalMatrix const &localMat)
Insert an element-matrix with row-indices == col-indices.
Definition: MatrixBackend.hpp:78
void zeroRows(std::vector< RowInd > const &rowLocalInd, bool diag)
Set all entries of the specified rows to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:114
void zeroRowsColumns(std::vector< RowInd > const &rowLocalInd, bool diag)
Set all entries of the specified rows and columns to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:152
void finish()
Finish assembly. Must be called before matrix can be used in a KSP.
Definition: MatrixBackend.hpp:238
BaseMatrix & matrix()
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:56
::Mat BaseMatrix
The matrix type of the underlying base matrix.
Definition: MatrixBackend.hpp:29
void zeroRowsColumns(std::vector< RowInd > const &rowLocalInd, bool diag, PetscVector< DofMap > const &x, PetscVector< DofMap > &b)
Set all entries of the specified rows and columns to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:163
void zeroRows(std::vector< RowInd > const &rowLocalInd, std::vector< ColInd > const &colLocalInd, bool diag)
Set all entries of the specified rows to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:126
void zeroRowsColumns(std::vector< RowInd > const &rowLocalInd, std::vector< ColInd > const &colLocalInd, bool diag)
Set all entries of the specified rows and columns to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:174
PetscMatrix(DofMap const &rowDofMap, DofMap const &colDofMap)
Constructor. Constructs new BaseMatrix.
Definition: MatrixBackend.hpp:39
void zeroRowsColumns(std::vector< RowInd > const &rowLocalInd, std::vector< ColInd > const &colLocalInd, bool diag, PetscVector< DofMap > const &x, PetscVector< DofMap > &b)
Set all entries of the specified rows and columns to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:181
PetscInt size_type
The index/size - type.
Definition: MatrixBackend.hpp:35
void scatter(RowLocalInd const &rowLocalInd, ColLocalInd const &colLocalInd, LocalMatrix const &localMat)
Insert an element-matrix.
Definition: MatrixBackend.hpp:96
PetscScalar value_type
The type of the elements of the DOFMatrix.
Definition: MatrixBackend.hpp:32
void init(MatrixNnzStructure const &pattern)
Create and initialize the matrix.
Definition: MatrixBackend.hpp:187
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:23
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:90
Definition: Constraints.hpp:14