AMDiS  0.1
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/linearalgebra/SymmetryStructure.hpp>
13 
14 namespace AMDiS
15 {
16  template <class> class Constraints;
17 
19  template <class DofMap>
21  {
22  template <class> friend class Constraints;
23 
24  public:
26  using BaseMatrix = ::Mat;
27 
29  using value_type = PetscScalar;
30 
32  using size_type = PetscInt;
33 
34  public:
36  template <class Basis>
37  explicit PetscMatrix(Basis const& rowBasis, Basis const& /*colBasis*/)
38  : dofMap_(rowBasis.communicator().dofMap())
39  {}
40 
41  // disable copy and move operations
42  PetscMatrix(PetscMatrix const&) = delete;
43  PetscMatrix(PetscMatrix&&) = delete;
44  PetscMatrix& operator=(PetscMatrix const&) = delete;
45  PetscMatrix& operator=(PetscMatrix&&) = delete;
46 
47  ~PetscMatrix()
48  {
49  destroy();
50  }
51 
54  {
55  return matrix_;
56  }
57 
59  BaseMatrix const& matrix() const
60  {
61  return matrix_;
62  }
63 
65  template <class RowIndex, class ColIndex>
66  void insert(RowIndex const& r, ColIndex const& c, PetscScalar value)
67  {
68  PetscInt row = dofMap_.global(r);
69  PetscInt col = dofMap_.global(c);
70  MatSetValue(matrix_, row, col, value, ADD_VALUES);
71  }
72 
74  template <class LocalInd, class LocalMatrix>
75  void scatter(LocalInd const& localInd, LocalMatrix const& localMat)
76  {
77  thread_local std::vector<PetscInt> idx;
78  idx.resize(localInd.size());
79 
80  // create a vector of global indices from the local indices using the local-to-global map
81  std::transform(localInd.begin(), localInd.end(), idx.begin(),
82  [this](auto const& mi) { return dofMap_.global(mi); });
83 
84  MatSetValues(matrix_, idx.size(), idx.data(), idx.size(), idx.data(), localMat.data(), ADD_VALUES);
85  }
86 
88  template <class RowLocalInd, class ColLocalInd, class LocalMatrix>
89  void scatter(RowLocalInd const& rowLocalInd, ColLocalInd const& colLocalInd, LocalMatrix const& localMat)
90  {
91  thread_local std::vector<PetscInt> ri;
92  thread_local std::vector<PetscInt> ci;
93  ri.resize(rowLocalInd.size());
94  ci.resize(colLocalInd.size());
95 
96  // create vectors of global indices from the local indices using the local-to-global map
97  std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
98  [this](auto const& mi) { return dofMap_.global(mi); });
99  std::transform(colLocalInd.begin(), colLocalInd.end(), ci.begin(),
100  [this](auto const& mi) { return dofMap_.global(mi); });
101 
102  MatSetValues(matrix_, ri.size(), ri.data(), ci.size(), ci.data(), localMat.data(), ADD_VALUES);
103  }
104 
106  template <class Pattern>
107  void init(Pattern const& pattern)
108  {
109  Dune::Timer t;
110 
111  // destroy an old matrix if created before
112  destroy();
113  info(3, " destroy old matrix needed {} seconds", t.elapsed());
114  t.reset();
115 
116  // create a MATAIJ or MATSEQAIJ matrix with provided sparsity pattern
117  MatCreateAIJ(comm(),
118  dofMap_.localSize(), dofMap_.localSize(),
119  dofMap_.globalSize(), dofMap_.globalSize(),
120  0, pattern.d_nnz().data(), 0, pattern.o_nnz().data(), &matrix_);
121 
122  // keep sparsity pattern even if we delete a row / column with e.g. MatZeroRows
123  MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
124 
125  // set symmetry properties of the matrix
126  switch (pattern.symmetry()) {
127  case SymmetryStructure::spd:
128  MatSetOption(matrix_, MAT_SPD, PETSC_TRUE);
129  break;
130  case SymmetryStructure::symmetric:
131  MatSetOption(matrix_, MAT_SYMMETRIC, PETSC_TRUE);
132  break;
133  case SymmetryStructure::hermitian:
134  MatSetOption(matrix_, MAT_HERMITIAN, PETSC_TRUE);
135  break;
136  case SymmetryStructure::structurally_symmetric:
137  MatSetOption(matrix_, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
138  break;
139  default:
140  /* do nothing */
141  break;
142  }
143 
144  info(3, " create new matrix needed {} seconds", t.elapsed());
145  t.reset();
146 
147  initialized_ = true;
148  }
149 
151  void init()
152  {
153  MatZeroEntries(matrix_);
154  initialized_ = true;
155  }
156 
158  void finish()
159  {
160  Dune::Timer t;
161  MatAssemblyBegin(matrix_, MAT_FINAL_ASSEMBLY);
162  MatAssemblyEnd(matrix_, MAT_FINAL_ASSEMBLY);
163  info(3, " finish matrix assembling needed {} seconds", t.elapsed());
164  }
165 
167  std::size_t nnz() const
168  {
169  MatInfo info;
170  MatGetInfo(matrix_, MAT_LOCAL, &info);
171  return std::size_t(info.nz_used);
172  }
173 
174  private:
175  // An MPI Communicator or PETSC_COMM_SELF
176  MPI_Comm comm() const
177  {
178 #if HAVE_MPI
179  return Environment::comm();
180 #else
181  return PETSC_COMM_SELF;
182 #endif
183  }
184 
185  // Destroy the matrix if initialized before.
186  void destroy()
187  {
188  if (initialized_)
189  MatDestroy(&matrix_);
190  }
191 
192  private:
193  // The local-to-global map
194  DofMap const& dofMap_;
195 
197  Mat matrix_;
198 
199  bool initialized_ = false;
200  };
201 
202 } // end namespace AMDiS
void finish()
Finish assembly. Must be called before matrix can be used in a KSP.
Definition: MatrixBackend.hpp:158
BaseMatrix const & matrix() const
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:59
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
void scatter(LocalInd const &localInd, LocalMatrix const &localMat)
Insert an element-matrix with row-indices == col-indices.
Definition: MatrixBackend.hpp:75
void init(Pattern const &pattern)
Create and initialize the matrix.
Definition: MatrixBackend.hpp:107
::Mat BaseMatrix
The matrix type of the underlying base matrix.
Definition: MatrixBackend.hpp:26
void scatter(RowLocalInd const &rowLocalInd, ColLocalInd const &colLocalInd, LocalMatrix const &localMat)
Insert an element-matrix.
Definition: MatrixBackend.hpp:89
The basic container that stores a base matrix.
Definition: MatrixBackend.hpp:20
PetscInt size_type
The index/size - type.
Definition: MatrixBackend.hpp:32
PetscScalar value_type
The type of the elements of the DOFMatrix.
Definition: MatrixBackend.hpp:29
void info(int noInfoLevel, std::string const &str, Args &&... args)
prints a message, if Environment::infoLevel() >= noInfoLevel
Definition: Output.hpp:105
PetscMatrix(Basis const &rowBasis, Basis const &)
Constructor. Constructs new BaseMatrix.
Definition: MatrixBackend.hpp:37
void init()
Reuse the matrix pattern and set all entries to zero.
Definition: MatrixBackend.hpp:151
BaseMatrix & matrix()
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:53
std::size_t nnz() const
Return the local number of nonzeros in the matrix.
Definition: MatrixBackend.hpp:167
Definition: Constraints.hpp:13
void insert(RowIndex const &r, ColIndex const &c, PetscScalar value)
Insert a single value into the matrix.
Definition: MatrixBackend.hpp:66
static Dune::MPIHelper::MPICommunicator comm()
Return the MPI_Comm object (or a fake communicator)
Definition: Environment.hpp:86