DOFMatrix.hpp 3.79 KB
Newer Older
1
2
#pragma once

3
#include <list>
4
#include <memory>
5
6
#include <string>
#include <vector>
7
8
9
10
11
12

#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#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
15
#include <amdis/linearalgebra/Common.hpp>
#include <amdis/linearalgebra/DOFMatrixBase.hpp>
16
17
18

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

28
    /// The type of the elements of the DOFMatrix
29
    using value_type = ValueType;
30

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

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

37
38
  public:
    /// Constructor. Constructs new BaseMatrix.
39
    MtlMatrix() = default;
40

41
    /// Return a reference to the data-matrix \ref matrix
42
    BaseMatrix& matrix()
43
    {
44
      assert( !inserter_ );
45
      return matrix_;
46
    }
47

48
    /// Return a reference to the data-matrix \ref matrix
49
    BaseMatrix const& matrix() const
50
    {
51
      assert( !inserter_ );
52
      return matrix_;
53
    }
54

55

56
    /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
57
    /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
58
    /// be created by \ref init and must be closed by \ref finish after insertion.
59
    void insert(size_type r, size_type c, value_type const& value)
60
    {
61
      test_exit_dbg(inserter_, "Inserter not initilized!");
62
      test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
63
          "Indices out of range [0,{})x[0,{})", num_rows(matrix_), num_cols(matrix_));
64
      (*inserter_)[r][c] += value;
65
66
    }

67

68
    /// Create inserter. Assumes that no inserter is currently active on this matrix.
69
70
71
    template <class RowBasis, class ColBasis>
    void init(RowBasis const& rowBasis, ColBasis const& colBasis,
              bool prepareForInsertion)
72
    {
73
      test_exit(!inserter_, "Matrix already in insertion mode!");
74

75
      calculateSlotSize();
76
      matrix_.change_dim(rowBasis.dimension(), colBasis.dimension());
77
      if (prepareForInsertion) {
78
79
        set_to_zero(matrix_);
        inserter_ = new Inserter(matrix_, slotSize_);
80
      }
81
82
    }

83
    /// Delete inserter -> finish insertion. Must be called in order to fill the
84
    /// final construction of the matrix.
85
86
    void finish()
    {
87
      delete inserter_;
88
      inserter_ = nullptr;
89
    }
90

91
    /// Return the number of nonzeros in the matrix
92
    std::size_t nnz() const
93
    {
94
      return matrix_.nnz();
95
    }
96

97
98
99
  protected:
    // Estimates the slot-size used in the inserter to reserve enough space per row.
    void calculateSlotSize()
100
    {
101
      slotSize_ = 0;
102

103
104
      if (num_rows(matrix_) != 0)
        slotSize_ = int(double(matrix_.nnz()) / num_rows(matrix_) * 1.2);
105
106
      if (slotSize_ < MIN_NNZ_PER_ROW)
        slotSize_ = MIN_NNZ_PER_ROW;
107
    }
108

109
  private:
110
    /// The minimal number of nonzeros per row
111
    static constexpr int MIN_NNZ_PER_ROW = 50;
112

113
    /// The data-matrix (can hold a new BaseMatrix or a pointer to external data
114
    BaseMatrix matrix_;
115

116
    /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
117
    Inserter* inserter_ = nullptr;
118

119
    /// The size of the slots used to insert values per row
120
    int slotSize_ = MIN_NNZ_PER_ROW;
121
122
  };

123
124
125
  template <class RowBasisType, class ColBasisType, class ValueType = double>
  using DOFMatrix = DOFMatrixBase<RowBasisType, ColBasisType, MtlMatrix<ValueType>>;

126
} // end namespace AMDiS