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