dune-istl  2.8.0
foreach.hh
Go to the documentation of this file.
1 #pragma once
2 
3 #include<type_traits>
4 #include<utility>
5 #include<cassert>
6 
7 #include<dune/common/diagonalmatrix.hh>
8 #include<dune/common/hybridutilities.hh>
9 #include<dune/common/indices.hh>
10 
13 
14 namespace Dune{
15 
16  namespace Impl {
17 
18  // stolen from dune-functions: we call a type "scalar" if it does not support index accessing
19  template<class C>
20  using StaticIndexAccessConcept = decltype(std::declval<C>()[Dune::Indices::_0]);
21 
22  template<class C>
23  using IsScalar = std::negation<Dune::Std::is_detected<StaticIndexAccessConcept, std::remove_reference_t<C>>>;
24 
25  // Type trait for matrix types that supports
26  // - iteration done row-wise
27  // - sparse iteration over nonzero entries
28  template <class T>
29  struct IsRowMajorSparse : std::false_type {};
30 
31  // This is supported by the following matrix types:
32  template <class B, class A>
33  struct IsRowMajorSparse<BCRSMatrix<B,A>> : std::true_type {};
34 
35  template <class K, int n>
36  struct IsRowMajorSparse<DiagonalMatrix<K,n>> : std::true_type {};
37 
38  template <class K, int n>
39  struct IsRowMajorSparse<ScaledIdentityMatrix<K,n>> : std::true_type {};
40 
41 
42  template <class Matrix>
43  auto rows(Matrix const& /*matrix*/, PriorityTag<2>) -> std::integral_constant<std::size_t, Matrix::N()> { return {}; }
44 
45  template <class Matrix>
46  auto cols(Matrix const& /*matrix*/, PriorityTag<2>) -> std::integral_constant<std::size_t, Matrix::M()> { return {}; }
47 
48  template <class Matrix>
49  auto rows(Matrix const& matrix, PriorityTag<1>) -> decltype(matrix.N()) { return matrix.N(); }
50 
51  template <class Matrix>
52  auto cols(Matrix const& matrix, PriorityTag<1>) -> decltype(matrix.M()) { return matrix.M(); }
53 
54  template <class Vector>
55  auto size(Vector const& /*vector*/, PriorityTag<2>) -> std::integral_constant<std::size_t, Vector::size()> { return {}; }
56 
57  template <class Vector>
58  auto size(Vector const& vector, PriorityTag<1>) -> decltype(vector.size()) { return vector.size(); }
59 
60 
61  } // end namespace Impl
62 
63 namespace ForEach{
64 
65  template <class Matrix>
66  auto rows(Matrix const& matrix) { return Impl::rows(matrix, PriorityTag<5>{}); }
67 
68  template <class Matrix>
69  auto cols(Matrix const& matrix) { return Impl::cols(matrix, PriorityTag<5>{}); }
70 
71  template <class Vector>
72  auto size(Vector const& vector) { return Impl::size(vector, PriorityTag<5>{}); }
73 
74 } // namespace ForEach
75 
76 
77 
78 
91 template <class Vector, class F>
92 std::size_t flatVectorForEach(Vector&& vector, F&& f, std::size_t offset = 0)
93 {
94  using V = std::decay_t<Vector>;
95  if constexpr( Impl::IsScalar<V>::value )
96  {
97  f(vector, offset);
98  return 1;
99  }
100  else
101  {
102  std::size_t idx = 0;
103  Hybrid::forEach(Dune::range(ForEach::size(vector)), [&](auto i) {
104  idx += flatVectorForEach(vector[i], f, offset + idx);
105  });
106  return idx;
107  }
108 }
109 
110 
128 template <class Matrix, class F>
129 std::pair<std::size_t,std::size_t> flatMatrixForEach(Matrix&& matrix, F&& f, std::size_t rowOffset = 0, std::size_t colOffset = 0)
130 {
131  using M = std::decay_t<Matrix>;
132  if constexpr ( Impl::IsScalar<M>::value )
133  {
134  f(matrix,rowOffset,colOffset);
135  return {1,1};
136  }
137  else
138  {
139  // if M supports the IsRowMajorSparse type trait: iterate just over the nonzero entries and
140  // and compute the flat row/col size directly
141  if constexpr ( Impl::IsRowMajorSparse<M>::value )
142  {
143  using Block = std::decay_t<decltype(matrix[0][0])>;
144 
145  // find an existing block or at least try to create one
146  auto block = [&]{
147  for (auto const& row : matrix)
148  for (auto const& entry : row)
149  return entry;
150  return Block{};
151  }();
152 
153  // compute the scalar size of the block
154  auto [blockRows, blockCols] = flatMatrixForEach(block, [](...){});
155 
156  // check whether we have valid sized blocks
157  assert( ( blockRows!=0 or blockCols!=0 ) and "the block size can't be zero");
158 
159  for ( auto rowIt = matrix.begin(); rowIt != matrix.end(); rowIt++ )
160  {
161  auto&& row = *rowIt;
162  auto rowIdx = rowIt.index();
163  for ( auto colIt = row.begin(); colIt != row.end(); colIt++ )
164  {
165  auto&& entry = *colIt;
166  auto colIdx = colIt.index();
167  auto [ dummyRows, dummyCols ] = flatMatrixForEach(entry, f, rowOffset + rowIdx*blockRows, colOffset + colIdx*blockCols);
168  assert( dummyRows == blockRows and dummyCols == blockCols and "we need the same size of each block in this matrix type");
169  }
170  }
171 
172  return { matrix.N()*blockRows, matrix.M()*blockCols };
173  }
174  // all other matrix types are accessed index-wise with dynamic flat row/col counting
175  else
176  {
177  std::size_t r = 0, c = 0;
178  std::size_t nRows, nCols;
179 
180  Hybrid::forEach(Dune::range(ForEach::rows(matrix)), [&](auto i) {
181  c = 0;
182  Hybrid::forEach(Dune::range(ForEach::cols(matrix)), [&](auto j) {
183  std::tie(nRows,nCols) = flatMatrixForEach(matrix[i][j], f, rowOffset + r, colOffset + c);
184  c += nCols;
185  });
186  r += nRows;
187  });
188  return {r,c};
189  }
190  }
191 }
192 
193 } // namespace Dune
Implementation of the BCRSMatrix class.
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Definition: allocator.hh:9
std::size_t flatVectorForEach(Vector &&vector, F &&f, std::size_t offset=0)
Traverse a blocked vector and call a functor at each scalar entry.
Definition: foreach.hh:92
std::pair< std::size_t, std::size_t > flatMatrixForEach(Matrix &&matrix, F &&f, std::size_t rowOffset=0, std::size_t colOffset=0)
Traverse a blocked matrix and call a functor at each scalar entry.
Definition: foreach.hh:129
auto rows(Matrix const &matrix)
Definition: foreach.hh:66
auto cols(Matrix const &matrix)
Definition: foreach.hh:69
auto size(Vector const &vector)
Definition: foreach.hh:72
A generic dynamic dense matrix.
Definition: matrix.hh:559
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
size_type N() const
Return the number of rows.
Definition: matrix.hh:693