Loading...
Searching...
No Matches
foreach.hh
Go to the documentation of this file.
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#pragma once
#include<type_traits>
#include<utility>
#include<cassert>
#include<dune/common/std/type_traits.hh>
#include<dune/common/diagonalmatrix.hh>
#include<dune/common/hybridutilities.hh>
#include<dune/common/indices.hh>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/scaledidmatrix.hh>
namespace Dune{
namespace Impl {
// stolen from dune-functions: we call a type "scalar" if it does not support index accessing
template<class C>
using StaticIndexAccessConcept = decltype(std::declval<C>()[Dune::Indices::_0]);
template<class C>
using IsScalar = std::negation<Dune::Std::is_detected<StaticIndexAccessConcept, std::remove_reference_t<C>>>;
// Type trait for matrix types that supports
// - iteration done row-wise
// - sparse iteration over nonzero entries
template <class T>
struct IsRowMajorSparse : std::false_type {};
// This is supported by the following matrix types:
template <class B, class A>
struct IsRowMajorSparse<BCRSMatrix<B,A>> : std::true_type {};
template <class K, int n>
struct IsRowMajorSparse<DiagonalMatrix<K,n>> : std::true_type {};
template <class K, int n>
struct IsRowMajorSparse<ScaledIdentityMatrix<K,n>> : std::true_type {};
template <class Matrix>
auto rows(Matrix const& /*matrix*/, PriorityTag<2>) -> std::integral_constant<std::size_t, Matrix::N()> { return {}; }
template <class Matrix>
auto cols(Matrix const& /*matrix*/, PriorityTag<2>) -> std::integral_constant<std::size_t, Matrix::M()> { return {}; }
template <class Matrix>
auto rows(Matrix const& matrix, PriorityTag<1>) -> decltype(matrix.N()) { return matrix.N(); }
template <class Matrix>
auto cols(Matrix const& matrix, PriorityTag<1>) -> decltype(matrix.M()) { return matrix.M(); }
template <class Vector>
auto size(Vector const& /*vector*/, PriorityTag<2>) -> std::integral_constant<std::size_t, Vector::size()> { return {}; }
template <class Vector>
auto size(Vector const& vector, PriorityTag<1>) -> decltype(vector.size()) { return vector.size(); }
} // end namespace Impl
namespace ForEach{
template <class Matrix>
auto rows(Matrix const& matrix) { return Impl::rows(matrix, PriorityTag<5>{}); }
template <class Matrix>
auto cols(Matrix const& matrix) { return Impl::cols(matrix, PriorityTag<5>{}); }
template <class Vector>
auto size(Vector const& vector) { return Impl::size(vector, PriorityTag<5>{}); }
} // namespace ForEach
template <class Vector, class F>
std::size_t flatVectorForEach(Vector&& vector, F&& f, std::size_t offset = 0)
{
using V = std::decay_t<Vector>;
if constexpr( Impl::IsScalar<V>::value )
{
f(vector, offset);
return 1;
}
else
{
std::size_t idx = 0;
Hybrid::forEach(Dune::range(ForEach::size(vector)), [&](auto i) {
idx += flatVectorForEach(vector[i], f, offset + idx);
});
return idx;
}
}
template <class Matrix, class F>
std::pair<std::size_t,std::size_t> flatMatrixForEach(Matrix&& matrix, F&& f, std::size_t rowOffset = 0, std::size_t colOffset = 0)
{
using M = std::decay_t<Matrix>;
if constexpr ( Impl::IsScalar<M>::value )
{
f(matrix,rowOffset,colOffset);
return {1,1};
}
else
{
// if M supports the IsRowMajorSparse type trait: iterate just over the nonzero entries and
// and compute the flat row/col size directly
if constexpr ( Impl::IsRowMajorSparse<M>::value )
{
using Block = std::decay_t<decltype(matrix[0][0])>;
// find an existing block or at least try to create one
auto block = [&]{
for (auto const& row : matrix)
for (auto const& entry : row)
return entry;
return Block{};
}();
// compute the scalar size of the block
auto [blockRows, blockCols] = flatMatrixForEach(block, [](...){});
// check whether we have valid sized blocks
assert( ( blockRows!=0 or blockCols!=0 ) and "the block size can't be zero");
for (auto rowIt = matrix.begin(); rowIt != matrix.end(); ++rowIt)
{
auto&& row = *rowIt;
auto rowIdx = rowIt.index();
for (auto colIt = row.begin(); colIt != row.end(); colIt++)
{
auto&& entry = *colIt;
auto colIdx = colIt.index();
#ifndef NDEBUG
// only instantiate return value in debug mode (for the assert)
auto [ numRows, numCols ] =
#endif
flatMatrixForEach(entry, f, rowOffset + rowIdx*blockRows, colOffset + colIdx*blockCols);
assert( numRows == blockRows and numCols == blockCols and "we need the same size of each block in this matrix type");
}
}
return { matrix.N()*blockRows, matrix.M()*blockCols };
}
// all other matrix types are accessed index-wise with dynamic flat row/col counting
else
{
std::size_t r = 0, c = 0;
std::size_t nRows, nCols;
Hybrid::forEach(Dune::range(ForEach::rows(matrix)), [&](auto i) {
c = 0;
Hybrid::forEach(Dune::range(ForEach::cols(matrix)), [&](auto j) {
std::tie(nRows,nCols) = flatMatrixForEach(matrix[i][j], f, rowOffset + r, colOffset + c);
c += nCols;
});
r += nRows;
});
return {r,c};
}
}
}
} // namespace Dune
bcrsmatrix.hh
Implementation of the BCRSMatrix class.
scaledidmatrix.hh
This file implements a quadratic matrix of fixed size which is a multiple of the identity.
Dune
Definition allocator.hh:11
Dune::flatMatrixForEach
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
Dune::flatVectorForEach
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
Dune::ForEach
Definition foreach.hh:66
Dune::ForEach::rows
auto rows(Matrix const &matrix)
Definition foreach.hh:69
Dune::ForEach::cols
auto cols(Matrix const &matrix)
Definition foreach.hh:72
Dune::ForEach::size
auto size(Vector const &vector)
Definition foreach.hh:75
Dune::Matrix
A generic dynamic dense matrix.
Definition matrix.hh:561
Dune::Matrix::M
size_type M() const
Return the number of columns.
Definition matrix.hh:696
Dune::Matrix::N
size_type N() const
Return the number of rows.
Definition matrix.hh:691