MatrixBackend.hpp 5.58 KB
Newer Older
1
2
3
4
5
6
7
8
#pragma once

#include <algorithm>
#include <memory>
#include <vector>

#include <petscmat.h>

Praetorius, Simon's avatar
Praetorius, Simon committed
9
10
#include <dune/common/timer.hh>

11
12
13
14
15
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>

namespace AMDiS
{
16
17
  template <class> class Constraints;

18
  /// \brief The basic container that stores a base matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
19
20
  template <class DofMap>
  class PetscMatrix
21
  {
22
    template <class> friend class Constraints;
23
24
25
26
27
28
29
30
31
32
33
34
35

  public:
    /// The matrix type of the underlying base matrix
    using BaseMatrix = ::Mat;

    /// The type of the elements of the DOFMatrix
    using value_type = PetscScalar;

    /// The index/size - type
    using size_type = PetscInt;

  public:
    /// Constructor. Constructs new BaseMatrix.
36
    template <class Basis>
Praetorius, Simon's avatar
Praetorius, Simon committed
37
38
    explicit PetscMatrix(Basis const& rowBasis, Basis const& /*colBasis*/)
      : dofMap_(rowBasis.communicator().dofMap())
39
40
41
    {}

    // disable copy and move operations
Praetorius, Simon's avatar
Praetorius, Simon committed
42
43
44
45
    PetscMatrix(PetscMatrix const&) = delete;
    PetscMatrix(PetscMatrix&&) = delete;
    PetscMatrix& operator=(PetscMatrix const&) = delete;
    PetscMatrix& operator=(PetscMatrix&&) = delete;
46

Praetorius, Simon's avatar
Praetorius, Simon committed
47
    ~PetscMatrix()
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
    {
      destroy();
    }

    /// Return a reference to the data-matrix \ref matrix
    BaseMatrix& matrix()
    {
      return matrix_;
    }

    /// Return a reference to the data-matrix \ref matrix
    BaseMatrix const& matrix() const
    {
      return matrix_;
    }

    /// Insert a single value into the matrix
    template <class RowIndex, class ColIndex>
    void insert(RowIndex const& r, ColIndex const& c, PetscScalar value)
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
68
69
      PetscInt row = dofMap_.global(r);
      PetscInt col = dofMap_.global(c);
70
71
72
73
74
75
76
77
78
79
80
81
      MatSetValue(matrix_, row, col, value, ADD_VALUES);
    }

    /// Insert an element-matrix with row-indices == col-indices
    template <class LocalInd, class LocalMatrix>
    void scatter(LocalInd const& localInd, LocalMatrix const& localMat)
    {
      thread_local std::vector<PetscInt> idx;
      idx.resize(localInd.size());

      // create a vector of global indices from the local indices using the local-to-global map
      std::transform(localInd.begin(), localInd.end(), idx.begin(),
Praetorius, Simon's avatar
Praetorius, Simon committed
82
        [this](auto const& mi) { return dofMap_.global(mi); });
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

      MatSetValues(matrix_, idx.size(), idx.data(), idx.size(), idx.data(), localMat.data(), ADD_VALUES);
    }

    /// Insert an element-matrix
    template <class RowLocalInd, class ColLocalInd, class LocalMatrix>
    void scatter(RowLocalInd const& rowLocalInd, ColLocalInd const& colLocalInd, LocalMatrix const& localMat)
    {
      thread_local std::vector<PetscInt> ri;
      thread_local std::vector<PetscInt> ci;
      ri.resize(rowLocalInd.size());
      ci.resize(colLocalInd.size());

      // create vectors of global indices from the local indices using the local-to-global map
      std::transform(rowLocalInd.begin(), rowLocalInd.end(), ri.begin(),
Praetorius, Simon's avatar
Praetorius, Simon committed
98
        [this](auto const& mi) { return dofMap_.global(mi); });
99
      std::transform(colLocalInd.begin(), colLocalInd.end(), ci.begin(),
Praetorius, Simon's avatar
Praetorius, Simon committed
100
        [this](auto const& mi) { return dofMap_.global(mi); });
101
102
103
104
105

      MatSetValues(matrix_, ri.size(), ri.data(), ci.size(), ci.data(), localMat.data(), ADD_VALUES);
    }

    /// Create and initialize the matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
106
    template <class Pattern>
107
    void init(Pattern const& pattern)
108
109
110
111
112
113
114
115
116
117
    {
      Dune::Timer t;

      // destroy an old matrix if created before
      destroy();
      info(3, "  destroy old matrix needed {} seconds", t.elapsed());
      t.reset();

      // create a MATAIJ or MATSEQAIJ matrix with provided sparsity pattern
      MatCreateAIJ(comm(),
Praetorius, Simon's avatar
Praetorius, Simon committed
118
119
120
        dofMap_.localSize(), dofMap_.localSize(),
        dofMap_.globalSize(), dofMap_.globalSize(),
        0, pattern.d_nnz().data(), 0, pattern.o_nnz().data(), &matrix_);
121
122
123
124

      // keep sparsity pattern even if we delete a row / column with e.g. MatZeroRows
      MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);

Praetorius, Simon's avatar
Praetorius, Simon committed
125
      // set symmetry properties of the matrix
126
      switch (pattern.symmetry()) {
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
        case SymmetryStructure::spd:
          MatSetOption(matrix_, MAT_SPD, PETSC_TRUE);
          break;
        case SymmetryStructure::symmetric:
          MatSetOption(matrix_, MAT_SYMMETRIC, PETSC_TRUE);
          break;
        case SymmetryStructure::hermitian:
          MatSetOption(matrix_, MAT_HERMITIAN, PETSC_TRUE);
          break;
        case SymmetryStructure::structurally_symmetric:
          MatSetOption(matrix_, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE);
          break;
        default:
          /* do nothing */
          break;
      }

      info(3, "  create new matrix needed {} seconds", t.elapsed());
      t.reset();

      initialized_ = true;
    }

    /// Finish assembly. Must be called before matrix can be used in a KSP
    void finish()
    {
      Dune::Timer t;
      MatAssemblyBegin(matrix_, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matrix_, MAT_FINAL_ASSEMBLY);
      info(3, "  finish matrix assembling needed {} seconds", t.elapsed());
    }

    /// Return the local number of nonzeros in the matrix
    std::size_t nnz() const
    {
      MatInfo info;
      MatGetInfo(matrix_, MAT_LOCAL, &info);
      return std::size_t(info.nz_used);
    }

  private:
    // An MPI Communicator or PETSC_COMM_SELF
    MPI_Comm comm() const
    {
#if HAVE_MPI
      return Environment::comm();
#else
      return PETSC_COMM_SELF;
#endif
    }

    // Destroy the matrix if initialized before.
    void destroy()
    {
      if (initialized_)
        MatDestroy(&matrix_);
    }

  private:
Praetorius, Simon's avatar
Praetorius, Simon committed
186
187
    // The local-to-global map
    DofMap const& dofMap_;
188
189
190
191
192
193
194
195

    /// The data-matrix
    Mat matrix_;

    bool initialized_ = false;
  };

} // end namespace AMDiS