9#include <dune/common/timer.hh>
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>
19 template <
class>
class Constraints;
22 template <
class DofMap>
40 : rowDofMap_(rowDofMap)
41 , colDofMap_(colDofMap)
68 template <
class RowIndex,
class ColIndex>
69 void insert(RowIndex
const& r, ColIndex
const& c, PetscScalar value)
71 PetscInt row = rowDofMap_.global(r);
72 PetscInt col = colDofMap_.global(c);
73 MatSetValue(matrix_, row, col, value, ADD_VALUES);
77 template <
class LocalInd,
class LocalMatrix>
78 void scatter(LocalInd
const& localInd, LocalMatrix
const& localMat)
80 if (&rowDofMap_ == &colDofMap_) {
81 thread_local std::vector<PetscInt> idx;
82 idx.resize(localInd.size());
85 std::transform(localInd.begin(), localInd.end(), idx.begin(),
86 [
this](
auto const& mi) { return rowDofMap_.global(mi); });
88 MatSetValues(matrix_, idx.size(), idx.data(), idx.size(), idx.data(), localMat.data(), ADD_VALUES);
90 scatter(localInd, localInd, localMat);
95 template <
class RowLocalInd,
class ColLocalInd,
class LocalMatrix>
96 void scatter(RowLocalInd
const& rowLocalInd, ColLocalInd
const& colLocalInd, LocalMatrix
const& localMat)
98 thread_local std::vector<PetscInt> ri;
99 thread_local std::vector<PetscInt> ci;
100 ri.resize(rowLocalInd.size());
101 ci.resize(colLocalInd.size());
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); });
109 MatSetValues(matrix_, ri.size(), ri.data(), ci.size(), ci.data(), localMat.data(), ADD_VALUES);
113 template <
class RowInd>
114 void zeroRows(std::vector<RowInd>
const& rowLocalInd,
bool diag)
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); });
121 MatZeroRows(matrix_, ri.size(), ri.data(), diag ? 1.0 : 0.0, PETSC_NULL, PETSC_NULL);
125 template <
class RowInd,
class ColInd>
126 void zeroRows(std::vector<RowInd>
const& rowLocalInd, std::vector<ColInd>
const& colLocalInd,
bool diag)
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);
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);
151 template <
class RowInd>
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);
162 template <
class RowInd>
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());
173 template <
class RowInd,
class ColInd>
174 void zeroRowsColumns(std::vector<RowInd>
const& rowLocalInd, std::vector<ColInd>
const& colLocalInd,
bool diag)
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.");
180 template <
class RowInd,
class ColInd>
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.");
193 info(3,
" destroy old matrix needed {} seconds", t.elapsed());
198 rowDofMap_.localSize(), colDofMap_.localSize(),
199 rowDofMap_.globalSize(), colDofMap_.globalSize(),
200 0, pattern.d_nnz().data(), 0, pattern.
o_nnz().data(), &matrix_);
203 MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
207 case SymmetryStructure::spd:
208 MatSetOption(matrix_, MAT_SPD, PETSC_TRUE);
210 case SymmetryStructure::symmetric:
211 MatSetOption(matrix_, MAT_SYMMETRIC, PETSC_TRUE);
213 case SymmetryStructure::hermitian:
214 MatSetOption(matrix_, MAT_HERMITIAN, PETSC_TRUE);
216 case SymmetryStructure::structurally_symmetric:
217 MatSetOption(matrix_, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
224 info(3,
" create new matrix needed {} seconds", t.elapsed());
233 MatZeroEntries(matrix_);
241 MatAssemblyBegin(matrix_, MAT_FINAL_ASSEMBLY);
242 MatAssemblyEnd(matrix_, MAT_FINAL_ASSEMBLY);
243 info(3,
" finish matrix assembling needed {} seconds", t.elapsed());
250 MatGetInfo(matrix_, MAT_LOCAL, &info);
251 return std::size_t(info.nz_used);
255 MPI_Comm comm()
const
257 return rowDofMap_.comm();
265 MatDestroy(&matrix_);
270 DofMap
const& rowDofMap_;
271 DofMap
const& colDofMap_;
276 bool initialized_ =
false;
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