DUNE-FEM (unstable)

foreach.hh
1// SPDX-FileCopyrightText: Copyright © 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>
11#include<dune/common/hybridutilities.hh>
12#include<dune/common/indices.hh>
13
16
17namespace 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
66namespace 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
94template <class Vector, class F>
95std::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
131template <class Matrix, class F>
132std::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#ifndef NDEBUG
171 // only instantiate return value in debug mode (for the assert)
172 auto [ numRows, numCols ] =
173#endif
174 flatMatrixForEach(entry, f, rowOffset + rowIdx*blockRows, colOffset + colIdx*blockCols);
175 assert( numRows == blockRows and numCols == blockCols and "we need the same size of each block in this matrix type");
176 }
177 }
178
179 return { matrix.N()*blockRows, matrix.M()*blockCols };
180 }
181 // all other matrix types are accessed index-wise with dynamic flat row/col counting
182 else
183 {
184 std::size_t r = 0, c = 0;
185 std::size_t nRows, nCols;
186
187 Hybrid::forEach(Dune::range(ForEach::rows(matrix)), [&](auto i) {
188 c = 0;
189 Hybrid::forEach(Dune::range(ForEach::cols(matrix)), [&](auto j) {
190 std::tie(nRows,nCols) = flatMatrixForEach(matrix[i][j], f, rowOffset + r, colOffset + c);
191 c += nCols;
192 });
193 r += nRows;
194 });
195 return {r,c};
196 }
197 }
198}
199
200} // namespace Dune
Implementation of the BCRSMatrix class.
A generic dynamic dense matrix.
Definition: matrix.hh:561
size_type M() const
Return the number of columns.
Definition: matrix.hh:696
size_type N() const
Return the number of rows.
Definition: matrix.hh:691
This file implements a quadratic diagonal matrix of fixed size.
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
Dune namespace.
Definition: alignedallocator.hh:13
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
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
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
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)