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

#include <list>
#include <memory>
#include <string>
#include <vector>

Praetorius, Simon's avatar
Praetorius, Simon committed
8
9
#include <Eigen/SparseCore>

10
11
12
13
14
15
16
17
#include <dune/common/timer.hh>

#include <amdis/Output.hpp>

namespace AMDiS
{
  /// \brief The basic container that stores a base matrix and a corresponding
  /// row/column feSpace.
Praetorius, Simon's avatar
Praetorius, Simon committed
18
19
  template <class T, int Orientation = Eigen::RowMajor>
  class EigenSparseMatrix
20
21
22
  {
  public:
    /// The matrix type of the underlying base matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
23
    using BaseMatrix = Eigen::SparseMatrix<T, Orientation>;
24
25

    /// The type of the elements of the DOFMatrix
26
    using value_type = typename BaseMatrix::Scalar;
27
28
29
30
31
32

    /// The index/size - type
    using size_type = typename BaseMatrix::Index;

  public:
    /// Constructor. Constructs new BaseMatrix.
33
    template <class Basis>
Praetorius, Simon's avatar
Praetorius, Simon committed
34
    explicit EigenSparseMatrix(Basis const&, Basis const&) {}
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

    /// 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_;
    }


    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
    /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
    /// be created by \ref init and must be closed by \ref finish after insertion.
    void insert(size_type r, size_type c, value_type const& value)
    {
      test_exit_dbg(r < matrix_.rows() && c < matrix_.cols(),
          "Indices out of range [0,{})x[0,{})", matrix_.rows(), matrix_.cols());
      tripletList_.emplace_back(r, c, value);
    }

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
    template <class Ind, class LocalMat>
    void scatter(Ind const& idx, LocalMat const& mat)
    {
      scatter(idx, idx, mat);
    }

    template <class RowInd, class ColInd, class LocalMat>
    void scatter(RowInd const& rows, ColInd const& cols, LocalMat const& mat)
    {
      //tripletList_.reserve(tripletList_.size() + rows.size()*cols.size());
      for (size_type i = 0; i < size_type(rows.size()); ++i)
        for (size_type j = 0; j < size_type(cols.size()); ++j)
          tripletList_.emplace_back(rows[i], cols[j], mat[i][j]);
    }

74
75

    /// Create inserter. Assumes that no inserter is currently active on this matrix.
Praetorius, Simon's avatar
Praetorius, Simon committed
76
    template <class Pattern>
77
    void init(Pattern const& pattern)
78
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
79
      matrix_.resize(pattern.rows(), pattern.cols());
80
      matrix_.setZero();
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
    }

    /// Delete inserter -> finish insertion. Must be called in order to fill the
    /// final construction of the matrix.
    void finish()
    {
      matrix_.setFromTriplets(tripletList_.begin(), tripletList_.end());
      matrix_.makeCompressed();

      tripletList_.clear(); // NOTE: this does not free the memory
    }

    /// Return the number of nonzeros in the matrix
    std::size_t nnz() const
    {
      return matrix_.nonZeros();
    }

  private:
    /// The data-matrix
    BaseMatrix matrix_;

    /// A list of row,col indices and values
104
    std::vector<Eigen::Triplet<value_type, Eigen::Index>> tripletList_;
105
106
107
  };

} // end namespace AMDiS