AMDiS  2.10
The Adaptive Multi-Dimensional Simulation Toolbox
MatrixBackend.hpp
1 #pragma once
2 
3 #include <list>
4 #include <memory>
5 #include <string>
6 #include <vector>
7 
8 #include <boost/numeric/mtl/matrix/compressed2D.hpp>
9 #include <boost/numeric/mtl/matrix/inserter.hpp>
10 #include <boost/numeric/mtl/utility/property_map.hpp>
11 #include <boost/numeric/mtl/utility/range_wrapper.hpp>
12 
13 #include <amdis/Output.hpp>
14 #include <amdis/linearalgebra/SymmetryStructure.hpp>
15 #include <amdis/linearalgebra/mtl/SlotSize.hpp>
16 #include <amdis/linearalgebra/mtl/VectorBackend.hpp>
17 
18 namespace AMDiS
19 {
20  class DefaultIndexDistribution;
21 
23  template <class T>
25  {
26  public:
28  using BaseMatrix = mtl::compressed2D<T>;
29 
31  using value_type = typename BaseMatrix::value_type;
32 
34  using size_type = typename BaseMatrix::size_type;
35 
36  private:
38  using Inserter = mtl::mat::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
39 
41  using Pattern = SlotSize;
42 
43  public:
46  {}
47 
50  {
51  assert( !inserter_ );
52  return matrix_;
53  }
54 
56  BaseMatrix const& matrix() const
57  {
58  assert( !inserter_ );
59  return matrix_;
60  }
61 
64  void init(Pattern const& pattern)
65  {
66  test_exit(!inserter_, "Matrix already in insertion mode!");
67 
68  std::size_t slotSize = nnz() > 0 ? 6*nnz() / (5*num_rows(matrix_))
69  : pattern.rowSizeEstimate();
70  matrix_.change_dim(pattern.rows(), pattern.cols());
71  set_to_zero(matrix_);
72 
73  symmetry_ = pattern.symmetry();
74  inserter_ = new Inserter(matrix_, slotSize);
75  }
76 
79  void init()
80  {
81  test_exit(!inserter_, "Matrix already in insertion mode!");
82 
83  std::size_t slotSize = 6*nnz() / (5*num_rows(matrix_));
84  set_to_zero(matrix_);
85 
86  inserter_ = new Inserter(matrix_, slotSize);
87  }
88 
89 
92  void finish()
93  {
94  delete inserter_;
95  inserter_ = nullptr;
96  }
97 
98 
102  void insert(size_type r, size_type c, value_type const& value)
103  {
104  test_exit_dbg(inserter_, "Inserter not initialized!");
105  test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
106  "Indices out of range [0,{})x[0,{})", num_rows(matrix_), num_cols(matrix_));
107  if (value != value_type(0) || r == c)
108  (*inserter_)[r][c] += value;
109  }
110 
111  template <class Ind, class LocalMat>
112  void scatter(Ind const& idx, LocalMat const& mat)
113  {
114  scatter(idx, idx, mat);
115  }
116 
117  template <class RowInd, class ColInd, class LocalMat>
118  void scatter(RowInd const& rows, ColInd const& cols, LocalMat const& mat)
119  {
120  test_exit_dbg(inserter_, "Inserter not initialized!");
121  for (size_type i = 0; i < size_type(rows.size()); ++i)
122  for (size_type j = 0; j < size_type(cols.size()); ++j)
123  if (mat[i][j] != value_type(0) || i == j)
124  (*inserter_)[rows[i]][cols[j]] += mat[i][j];
125  }
126 
127  template <class RowInd, class ColInd>
128  void setUnitDiagonals(RowInd const& rowInd, ColInd const& colInd)
129  {
130  bool all_diagonals_found = true;
131  auto const& indexer = matrix_.indexer;
132  auto rowIt = rowInd.begin();
133  auto rowEnd = rowInd.end();
134  auto colIt = colInd.begin();
135  for (; rowIt != rowEnd; ++rowIt, ++colIt) {
136  auto pos = indexer(matrix_, *rowIt, *colIt);
137  if (pos)
138  matrix_.value_from_offset(pos.value()) = 1.0;
139  else {
140  all_diagonals_found = false;
141  break;
142  }
143  }
144 
145  if (!all_diagonals_found) {
146  std::size_t slotSize = 6*nnz() / (5*num_rows(matrix_));
147  mtl::mat::inserter<BaseMatrix, mtl::update_store<value_type> > ins(matrix_, slotSize);
148  auto rowIt = rowInd.begin();
149  auto rowEnd = rowInd.end();
150  auto colIt = colInd.begin();
151  for (; rowIt != rowEnd; ++rowIt, ++colIt)
152  ins(*rowIt,*colIt) << 1.0;
153  }
154  }
155 
157  template <class RowInd, class ColInd>
158  void zeroRows(RowInd const& rowInd, ColInd const& colInd, bool diag)
159  {
160  for (auto i : rowInd) {
161  size_type offset_begin = matrix_.ref_major()[i];
162  size_type offset_end = matrix_.ref_major()[i+1];
163  std::fill(matrix_.address_data()+offset_begin, matrix_.address_data()+offset_end, value_type(0));
164  }
165 
166  if (diag)
167  setUnitDiagonals(rowInd,colInd);
168  }
169 
170  template <class RowInd>
171  void zeroRows(RowInd const& rowInd, bool diag)
172  {
173  zeroRows(rowInd,rowInd,diag);
174  }
175 
177  template <class RowInd, class ColInd>
178  void zeroRowsColumnsImpl(RowInd const& rowInd, ColInd const& colInd, bool diag, MTLVector<T> const* x = nullptr, MTLVector<T>* b = nullptr)
179  {
180  assert(rowInd.size() == colInd.size());
181 
182  std::vector<bool> colBitVec(num_cols(matrix_), false);
183  for (std::size_t i = 0; i < colInd.size(); ++i)
184  colBitVec[colInd[i]] = true;
185 
186  zeroRows(rowInd, false);
187 
188  // Define the property maps
189  auto row = mtl::mat::row_map(matrix_);
190  auto col = mtl::mat::col_map(matrix_);
191  auto value = mtl::mat::value_map(matrix_);
192 
193  // iterate over the matrix
194  std::size_t slotSize = 0;
195  for (auto r : mtl::rows_of(matrix_)) { // rows or columns
196  std::size_t rowSize = 0;
197  for (auto i : mtl::nz_of(r)) { // non-zeros within
198  ++rowSize;
199  if (colBitVec[col(i)]) {
200  if (diag && x && b)
201  b->vector()[row(i)] -= value(i) * x->vector()[col(i)];
202  value(i, 0);
203  }
204  }
205  slotSize = std::max(slotSize, rowSize);
206  }
207 
208  // set diagonal entry
209  if (diag)
210  setUnitDiagonals(rowInd,colInd);
211  }
212 
213  template <class RowInd>
214  void zeroRowsColumns(RowInd const& rowInd, bool diag)
215  {
216  zeroRowsColumnsImpl(rowInd,rowInd,diag,nullptr,nullptr);
217  }
218 
219  template <class RowInd>
220  void zeroRowsColumns(RowInd const& rowInd, bool diag, MTLVector<T> const& x, MTLVector<T>& b)
221  {
222  zeroRowsColumnsImpl(rowInd,rowInd,diag,&x,&b);
223  }
224 
225  template <class RowInd, class ColInd>
226  void zeroRowsColumns(RowInd const& rowInd, ColInd const& colInd, bool diag)
227  {
228  zeroRowsColumnsImpl(rowInd,colInd,diag,nullptr,nullptr);
229  }
230 
231  template <class RowInd, class ColInd>
232  void zeroRowsColumns(RowInd const& rowInd, ColInd const& colInd, bool diag, MTLVector<T> const& x, MTLVector<T>& b)
233  {
234  zeroRowsColumnsImpl(rowInd,colInd,diag,&x,&b);
235  }
236 
238  std::size_t nnz() const
239  {
240  return matrix_.nnz();
241  }
242 
244  SymmetryStructure symmetry() const
245  {
246  return symmetry_;
247  }
248 
249  private:
251  BaseMatrix matrix_;
252 
254  Inserter* inserter_ = nullptr;
255 
257  SymmetryStructure symmetry_ = SymmetryStructure::unknown;
258  };
259 
260 } // end namespace AMDiS
BaseMatrix & matrix()
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:49
void init()
Definition: MatrixBackend.hpp:79
SymmetryStructure symmetry() const
Symmetry of the matrix entries.
Definition: MatrixBackend.hpp:244
std::size_t rowSizeEstimate() const
Estimate of the non-zeros per row.
Definition: SlotSize.hpp:41
mtl::compressed2D< T > BaseMatrix
The matrix type of the underlying base matrix.
Definition: MatrixBackend.hpp:28
BaseMatrix const & matrix() const
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:56
Definition: AdaptBase.hpp:6
std::size_t cols() const
Number of columns in the matrix.
Definition: SlotSize.hpp:35
typename BaseMatrix::size_type size_type
The index/size - type.
Definition: MatrixBackend.hpp:34
std::size_t rows() const
Number of rows in the matrix.
Definition: SlotSize.hpp:29
SymmetryStructure symmetry() const
Symmetry of the matrix entries.
Definition: SlotSize.hpp:47
void zeroRowsColumnsImpl(RowInd const &rowInd, ColInd const &colInd, bool diag, MTLVector< T > const *x=nullptr, MTLVector< T > *b=nullptr)
Set all entries of the specified rows and columns to zero and the diagonal element to diag ...
Definition: MatrixBackend.hpp:178
void init(Pattern const &pattern)
Definition: MatrixBackend.hpp:64
The basic container that stores a base matrix.
Definition: MatrixBackend.hpp:24
std::size_t nnz() const
Return the number of nonzeros in the matrix.
Definition: MatrixBackend.hpp:238
Definition: SlotSize.hpp:15
typename BaseMatrix::value_type value_type
The type of the elements of the DOFMatrix.
Definition: MatrixBackend.hpp:31
void insert(size_type r, size_type c, value_type const &value)
Returns an update-proxy of the inserter, to insert/update a value at position (r, c) in the matrix...
Definition: MatrixBackend.hpp:102
Dummy implementation for sequential index "distribution".
Definition: IndexDistribution.hpp:6
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:158
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:18
MTLSparseMatrix(DefaultIndexDistribution const &, DefaultIndexDistribution const &)
Constructor. Constructs new BaseMatrix.
Definition: MatrixBackend.hpp:45
void finish()
Definition: MatrixBackend.hpp:92