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

Praetorius, Simon's avatar
Praetorius, Simon committed
3
#include <list>
4
5
6
#include <string>
#include <memory>

Praetorius, Simon's avatar
Praetorius, Simon committed
7
#include <dune/istl/bcrsmatrix.hh>
8
9
#include <dune/istl/matrixindexset.hh>

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

namespace AMDiS
{
Praetorius, Simon's avatar
Praetorius, Simon committed
15
16
17
18
19
20
  template <class T, class = void>
  struct BlockMatrixType
  {
    using type = Dune::FieldMatrix<T,1,1>;
  };

21
  template <class T>
Praetorius, Simon's avatar
Praetorius, Simon committed
22
  struct BlockMatrixType<T, typename T::field_type>
23
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
24
25
    using type = T;
  };
26

Praetorius, Simon's avatar
Praetorius, Simon committed
27
  template <class T, class C>
Praetorius, Simon's avatar
Praetorius, Simon committed
28
29
30
  class ISTLBCRSMatrix
  {
  public:
31
    /// The matrix type of the underlying base matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
32
    using BaseMatrix = Dune::BCRSMatrix<typename BlockMatrixType<T>::type>;
33

Praetorius, Simon's avatar
Praetorius, Simon committed
34
35
36
    /// Communication type
    using Comm = C;

37
38
    /// The type of the elements of the DOFMatrix
    using value_type = typename BaseMatrix::block_type;
39

40
41
42
    /// The index/size - type
    using size_type = typename BaseMatrix::size_type;

43
  public:
44
    /// Constructor. Constructs new BaseVector.
45
    template <class Basis>
Praetorius, Simon's avatar
Praetorius, Simon committed
46
47
48
    explicit ISTLBCRSMatrix(Basis const& rowBasis, Basis const&)
      : comm_(&rowBasis.communicator())
    {}
49

50
    /// Return the data-vector \ref vector
51
    BaseMatrix const& matrix() const
52
    {
53
      return matrix_;
54
    }
55

56
    /// Return the data-vector \ref vector
57
    BaseMatrix& matrix()
58
    {
59
      return matrix_;
60
    }
61

Praetorius, Simon's avatar
Praetorius, Simon committed
62
63
    Comm const& comm() const { return *comm_; }

64
    /// create occupation pattern and apply it to the matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
65
    template <class Pattern>
66
    void init(Pattern const& pattern)
67
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
68
      pattern.applyTo(matrix_);
69
      initialized_ = true;
70
71
72
73
    }

    void finish()
    {
74
      initialized_ = false;
75
    }
76

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

    /// Insert a single value into the matrix (add to existing value)
    void insert(size_type r, size_type c, value_type const& value)
    {
      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
      test_exit_dbg( r < matrix_.N() && c < matrix_.M() ,
          "Indices out of range [0,{})x[0,{})", matrix_.N(), matrix_.M() );
      matrix_[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( initialized_, "Occupation pattern not initialized!");
      for (size_type i = 0; i < size_type(rows.size()); ++i)
        for (size_type j = 0; j < size_type(cols.size()); ++j)
          matrix_[rows[i]][cols[j]] += mat[i][j];
    }

102
    std::size_t nnz() const
103
    {
104
      return matrix_.nonzeroes();
105
    }
106

107
  private:
108
    BaseMatrix matrix_;
Praetorius, Simon's avatar
Praetorius, Simon committed
109
    Comm const* comm_;
110

111
    bool initialized_ = false;
112
113
114
  };

} // end namespace AMDiS