AMDiS  2.10
The Adaptive Multi-Dimensional Simulation Toolbox
MatrixBackend.hpp
1 #pragma once
2 
3 #include <list>
4 #include <string>
5 #include <memory>
6 
7 #include <dune/istl/bcrsmatrix.hh>
8 #include <dune/istl/matrixindexset.hh>
9 
10 #include <amdis/Output.hpp>
11 #include <amdis/linearalgebra/SparsityPattern.hpp>
12 #include <amdis/linearalgebra/SymmetryStructure.hpp>
13 #include <amdis/linearalgebra/istl/VectorBackend.hpp>
14 
15 namespace AMDiS
16 {
17  template <class T, class = void>
19  {
20  using type = Dune::FieldMatrix<T,1,1>;
21  };
22 
23  template <class T>
24  struct BlockMatrixType<T, typename T::field_type>
25  {
26  using type = T;
27  };
28 
29  template <class T, class CommType>
31  {
32  public:
34  using BaseMatrix = Dune::BCRSMatrix<typename BlockMatrixType<T>::type>;
35 
37  using Comm = CommType;
38 
40  using value_type = typename BaseMatrix::block_type;
41 
43  using size_type = typename BaseMatrix::size_type;
44 
45  public:
47  ISTLBCRSMatrix(Comm const& comm, Comm const&)
48  : comm_(comm)
49  {}
50 
52  BaseMatrix const& matrix() const
53  {
54  return matrix_;
55  }
56 
59  {
60  return matrix_;
61  }
62 
63  Comm const& comm() const { return comm_; }
64 
66  void init(SparsityPattern const& pattern)
67  {
68  pattern.applyTo(matrix_);
69  initialized_ = true;
70  }
71 
73  void init()
74  {
75  matrix_ = value_type(0);
76  initialized_ = true;
77  }
78 
79  void finish()
80  {
81  initialized_ = false;
82  }
83 
84 
86  void insert(size_type r, size_type c, value_type const& value)
87  {
88  test_exit_dbg( initialized_, "Occupation pattern not initialized!");
89  test_exit_dbg( r < matrix_.N() && c < matrix_.M() ,
90  "Indices out of range [0,{})x[0,{})", matrix_.N(), matrix_.M() );
91  matrix_[r][c] += value;
92  }
93 
94  template <class Ind, class LocalMat>
95  void scatter(Ind const& idx, LocalMat const& mat)
96  {
97  scatter(idx, idx, mat);
98  }
99 
100  template <class RowInd, class ColInd, class LocalMat>
101  void scatter(RowInd const& rows, ColInd const& cols, LocalMat const& mat)
102  {
103  test_exit_dbg( initialized_, "Occupation pattern not initialized!");
104  for (size_type i = 0; i < size_type(rows.size()); ++i)
105  for (size_type j = 0; j < size_type(cols.size()); ++j)
106  matrix_[rows[i]][cols[j]] += mat[i][j];
107  }
108 
110  template <class RowInd>
111  void zeroRows(RowInd const& rowInd, bool diag)
112  {
113  // loop over the matrix rows
114  for (size_type i : rowInd) {
115  auto cIt = matrix_[i].begin();
116  auto cEndIt = matrix_[i].end();
117  // loop over nonzero matrix entries in current row
118  for (; cIt != cEndIt; ++cIt)
119  *cIt = (diag && i == cIt.index() ? T(1) : T(0));
120  }
121  }
122 
124  template <class RowInd, class ColInd>
125  void zeroRows(RowInd const& rowInd, ColInd const& colInd, bool diag)
126  {
127  // loop over the matrix rows
128  auto rowIt = rowInd.begin();
129  auto rowEnd = rowInd.end();
130  auto colIt = colInd.begin();
131  for (; rowIt != rowEnd; ++rowIt,++colIt) {
132  auto cIt = matrix_[*rowIt].begin();
133  auto cEndIt = matrix_[*rowIt].end();
134  // loop over nonzero matrix entries in current row
135  for (; cIt != cEndIt; ++cIt)
136  *cIt = (diag && *colIt == cIt.index() ? T(1) : T(0));
137  }
138  }
139 
141  template <class RowInd, class ColInd>
142  void zeroRowsColumnsImpl(RowInd const& rowInd, ColInd const& colInd, bool diag, ISTLBlockVector<T> const* x = nullptr, ISTLBlockVector<T>* b = nullptr)
143  {
144  assert(rowInd.size() == colInd.size());
145 
146  std::vector<bool> colBitVec(matrix_.M(), false);
147  for (std::size_t i = 0; i < colInd.size(); ++i)
148  colBitVec[colInd[i]] = true;
149 
150  // loop over nonzero matrix entries in current row
151  for (size_type i : rowInd) {
152  auto cIt = matrix_[i].begin();
153  auto cEndIt = matrix_[i].end();
154  for (; cIt != cEndIt; ++cIt)
155  *cIt = T(0);
156  }
157 
158  // iterate over the matrix
159  for (size_type i = 0; i < matrix_.N(); ++i) {
160  auto cIt = matrix_[i].begin();
161  auto cEndIt = matrix_[i].end();
162  // loop over nonzero matrix entries in current row
163  for (; cIt != cEndIt; ++cIt) {
164  if (colBitVec[cIt.index()]) {
165  if (diag && x && b)
166  b->vector()[i] -= (*cIt) * x->vector()[cIt.index()];
167  *cIt = 0;
168  }
169  }
170  }
171 
172  // set diagonal entry
173  if (diag) {
174  for (std::size_t i = 0; i < rowInd.size(); ++i) {
175  test_exit_dbg(matrix_.exists(rowInd[i], colInd[i]),
176  "The entry (",rowInd[i],",",colInd[i],") does not yet exist. Must be inserted in the pattern.");
177  matrix_[rowInd[i]][colInd[i]] = 1;
178  }
179  }
180  }
181  template <class RowInd>
182  void zeroRowsColumns(RowInd const& rowInd, bool diag)
183  {
184  zeroRowsColumnsImpl(rowInd,rowInd,diag,nullptr,nullptr);
185  }
186 
187  template <class RowInd>
188  void zeroRowsColumns(RowInd const& rowInd, bool diag, ISTLBlockVector<T> const& x, ISTLBlockVector<T>& b)
189  {
190  zeroRowsColumnsImpl(rowInd,rowInd,diag,&x,&b);
191  }
192 
193  template <class RowInd, class ColInd>
194  void zeroRowsColumns(RowInd const& rowInd, ColInd const& colInd, bool diag)
195  {
196  zeroRowsColumnsImpl(rowInd,colInd,diag,nullptr,nullptr);
197  }
198 
199  template <class RowInd, class ColInd>
200  void zeroRowsColumns(RowInd const& rowInd, ColInd const& colInd, bool diag, ISTLBlockVector<T> const& x, ISTLBlockVector<T>& b)
201  {
202  zeroRowsColumnsImpl(rowInd,colInd,diag,&x,&b);
203  }
204 
206  std::size_t nnz() const
207  {
208  return matrix_.nonzeroes();
209  }
210 
211  private:
212  BaseMatrix matrix_;
213  Comm const& comm_;
214 
215  bool initialized_ = false;
216  };
217 
218 } // end namespace AMDiS
BaseMatrix const & matrix() const
Return the data-vector vector.
Definition: MatrixBackend.hpp:52
void insert(size_type r, size_type c, value_type const &value)
Insert a single value into the matrix (add to existing value)
Definition: MatrixBackend.hpp:86
void zeroRows(RowInd const &rowInd, bool diag)
Set all entries of the specified rows to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:111
Definition: AdaptBase.hpp:6
void init(SparsityPattern const &pattern)
Create occupation pattern and apply it to the matrix.
Definition: MatrixBackend.hpp:66
Dune::BCRSMatrix< typename BlockMatrixType< T >::type > BaseMatrix
The matrix type of the underlying base matrix.
Definition: MatrixBackend.hpp:34
CommType Comm
Communication type.
Definition: MatrixBackend.hpp:37
void zeroRows(RowInd const &rowInd, ColInd const &colInd, bool diag)
Set all entries of the specified rows to zero and the diagonal element to diag
Definition: MatrixBackend.hpp:125
ISTLBCRSMatrix(Comm const &comm, Comm const &)
Constructor. Constructs new BaseVector.
Definition: MatrixBackend.hpp:47
Definition: MatrixBackend.hpp:30
void applyTo(Matrix &matrix) const
Initialize a matrix with the non-zero structure.
Definition: SparsityPattern.hpp:75
typename BaseMatrix::size_type size_type
The index/size - type.
Definition: MatrixBackend.hpp:43
A general sparsity pattern implementation using the full pattern of the basis by adding all local ind...
Definition: SparsityPattern.hpp:14
std::size_t nnz() const
Return the number of entries allocated in the sparsity pattern of the matrix.
Definition: MatrixBackend.hpp:206
BaseMatrix & matrix()
Return the data-vector vector.
Definition: MatrixBackend.hpp:58
void init()
Set all entries to zero while keeping the occupation pattern intact.
Definition: MatrixBackend.hpp:73
void zeroRowsColumnsImpl(RowInd const &rowInd, ColInd const &colInd, bool diag, ISTLBlockVector< T > const *x=nullptr, ISTLBlockVector< T > *b=nullptr)
Set all entries of the specified rows and columns to zero and the diagonal element to diag ...
Definition: MatrixBackend.hpp:142
Definition: MatrixBackend.hpp:18
typename BaseMatrix::block_type value_type
The type of the elements of the DOFMatrix.
Definition: MatrixBackend.hpp:40
Definition: VectorBackend.hpp:26