MatrixBackend.hpp 3.98 KB
Newer Older
1
2
#pragma once

3
#include <list>
4
#include <memory>
5
6
#include <string>
#include <vector>
7

Praetorius, Simon's avatar
Praetorius, Simon committed
8
#include <boost/numeric/mtl/matrix/compressed2D.hpp>
9
10
11
12
#include <boost/numeric/mtl/matrix/inserter.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>

13
#include <amdis/Output.hpp>
14
#include <amdis/linearalgebra/SymmetryStructure.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
15
#include <amdis/linearalgebra/mtl/SlotSize.hpp>
16
17
18

namespace AMDiS
{
19
  /// \brief The basic container that stores a base matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
20
21
  template <class T>
  class MTLSparseMatrix
22
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
23
  public:
24
    /// The matrix type of the underlying base matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
25
    using BaseMatrix = mtl::compressed2D<T>;
26

27
    /// The type of the elements of the DOFMatrix
28
    using value_type = typename BaseMatrix::value_type;
29

30
31
32
    /// The index/size - type
    using size_type = typename BaseMatrix::size_type;

Praetorius, Simon's avatar
Praetorius, Simon committed
33
  private:
34
    /// The type of the inserter used when filling the matrix. NOTE: need not be public
35
    using Inserter = mtl::mat::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
36

Praetorius, Simon's avatar
Praetorius, Simon committed
37
38
39
    /// Type of the matrix non-zero pattern
    using Pattern = SlotSize;

40
41
  public:
    /// Constructor. Constructs new BaseMatrix.
42
    template <class Basis>
Praetorius, Simon's avatar
Praetorius, Simon committed
43
    explicit MTLSparseMatrix(Basis const& /*rowBasis*/, Basis const& /*colBasis*/) {}
44

45
    /// Return a reference to the data-matrix \ref matrix
46
    BaseMatrix& matrix()
47
    {
48
      assert( !inserter_ );
49
      return matrix_;
50
    }
51

52
    /// Return a reference to the data-matrix \ref matrix
53
    BaseMatrix const& matrix() const
54
    {
55
      assert( !inserter_ );
56
      return matrix_;
57
    }
58

59
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
Praetorius, Simon's avatar
Praetorius, Simon committed
60
    void init(Pattern const& pattern, SymmetryStructure symmetry)
61
    {
62
      test_exit(!inserter_, "Matrix already in insertion mode!");
63

Praetorius, Simon's avatar
Praetorius, Simon committed
64
65
      std::size_t slotSize = pattern.avgRowSize(matrix_);
      matrix_.change_dim(pattern.rows(), pattern.cols());
66
67
      set_to_zero(matrix_);

Praetorius, Simon's avatar
Praetorius, Simon committed
68
      inserter_ = new Inserter(matrix_, slotSize);
69
      symmetry_ = symmetry;
70
71
    }

72
    /// Delete inserter -> finish insertion. Must be called in order to fill the
73
    /// final construction of the matrix.
74
75
    void finish()
    {
76
      delete inserter_;
77
      inserter_ = nullptr;
78
    }
79

80

Praetorius, Simon's avatar
Praetorius, Simon committed
81
    /// \brief Returns an update-proxy of the inserter, to insert/update a value at
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    /// 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(inserter_, "Inserter not initilized!");
      test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
          "Indices out of range [0,{})x[0,{})", num_rows(matrix_), num_cols(matrix_));
      if (value != value_type(0) || r == c)
        (*inserter_)[r][c] += value;
    }

    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)
    {
      test_exit_dbg(inserter_, "Inserter not initilized!");
      for (size_type i = 0; i < size_type(rows.size()); ++i)
        for (size_type j = 0; j < size_type(cols.size()); ++j)
          if (mat[i][j] != value_type(0) || i == j)
            (*inserter_)[rows[i]][cols[j]] += mat[i][j];
    }

109
    /// Return the number of nonzeros in the matrix
110
    std::size_t nnz() const
111
    {
112
      return matrix_.nnz();
113
    }
114

Praetorius, Simon's avatar
Praetorius, Simon committed
115
    /// Symmetry of the matrix entries
116
117
118
119
120
    SymmetryStructure symmetry() const
    {
      return symmetry_;
    }

121
  private:
122
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
123
    BaseMatrix matrix_;
124

125
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
126
    Inserter* inserter_ = nullptr;
127

Praetorius, Simon's avatar
Praetorius, Simon committed
128
    /// Symmetry of the matrix entries
129
130
    SymmetryStructure symmetry_ = SymmetryStructure::unknown;
  };
131

132
} // end namespace AMDiS