Loading...
Searching...
No Matches
transpose.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_COMMON_TRANSPOSE_HH
#define DUNE_COMMON_TRANSPOSE_HH
#include <cstddef>
#include <functional>
#include <dune/common/boundschecking.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/common/referencehelper.hh>
#include <dune/common/dynmatrix.hh>
#include <dune/common/matrixconcepts.hh>
namespace Dune {
namespace Impl {
template<class M, bool = IsStaticSizeMatrix<M>::value>
struct TransposedDenseMatrixTraits
{
using type = Dune::FieldMatrix<typename FieldTraits<M>::field_type, M::cols, M::rows>;
};
template<class M>
struct TransposedDenseMatrixTraits<M, false>
{
using type = Dune::DynamicMatrix<typename FieldTraits<M>::field_type>;
};
template<class M>
using TransposedDenseMatrixTraits_t = typename TransposedDenseMatrixTraits<M>::type;
// CRTP mixin class to provide the static part of the interface,
// namely enums rows and cols.
template<class WM, bool staticSize = IsStaticSizeMatrix<WM>::value>
class TransposedMatrixWrapperMixin {};
template<class WM>
class TransposedMatrixWrapperMixin<WM, true>
{
public:
constexpr static int rows = WM::cols;
constexpr static int cols = WM::rows;
};
template<class M>
class TransposedMatrixWrapper;
} // namespace Impl
// Specialization of FieldTraits needs to be in namespace Dune::
template<class M>
struct FieldTraits< Impl::TransposedMatrixWrapper<M> >
{
using field_type = typename FieldTraits<ResolveRef_t<M>>::field_type;
using real_type = typename FieldTraits<ResolveRef_t<M>>::real_type;
};
namespace Impl {
// Specialize the traits `IsDenseMatrix` as true only if the wrapped matrix
// is a dense matrix.
template<class M>
struct IsDenseMatrix< TransposedMatrixWrapper<M> > :
public IsDenseMatrix<ResolveRef_t<M>> {};
// Wrapper representing the transposed of a matrix.
// Creating the wrapper does not compute anything
// but only serves for tagging the wrapped matrix
// for transposition. This class will store the matrix
// of type `Matrix` by value. To support reference-semantic,
// it supports using Matrix=std::reference_wrapper<OriginalMatrixType>.
template<class Matrix>
class TransposedMatrixWrapper :
public TransposedMatrixWrapperMixin<ResolveRef_t<Matrix>>
{
constexpr static bool hasStaticSize = IsStaticSizeMatrix<ResolveRef_t<Matrix>>::value;
public:
using WrappedMatrix = ResolveRef_t<Matrix>;
const WrappedMatrix& wrappedMatrix() const {
return resolveRef(matrix_);
}
using value_type = typename WrappedMatrix::value_type;
using field_type = typename FieldTraits<WrappedMatrix>::field_type;
using size_type = typename WrappedMatrix::size_type;
TransposedMatrixWrapper(Matrix&& matrix) : matrix_(std::move(matrix)) {}
TransposedMatrixWrapper(const Matrix& matrix) : matrix_(matrix) {}
private:
// helper proxy type to access elements by chained `operator[][]`
struct AccessProxy
{
const WrappedMatrix& wrappedMatrix_;
size_type i_;
decltype(auto) operator[] (size_type j) const
{
DUNE_ASSERT_BOUNDS(j < wrappedMatrix_.N());
return wrappedMatrix_[j][i_];
}
};
public:
auto operator[] (size_type i) const requires(Impl::IsDenseMatrix_v<WrappedMatrix>)
{
DUNE_ASSERT_BOUNDS(i < wrappedMatrix().M());
return AccessProxy{wrappedMatrix(),i};
}
template<class MatrixA,
std::enable_if_t<
((not Impl::IsFieldMatrix<MatrixA>::value) or (not hasStaticSize))
and Impl::IsDenseMatrix_v<MatrixA>, int> = 0>
friend auto operator* (const MatrixA& matrixA, const TransposedMatrixWrapper& matrixB)
{
using FieldA = typename FieldTraits<MatrixA>::field_type;
using FieldB = typename FieldTraits<TransposedMatrixWrapper>::field_type;
using Field = typename PromotionTraits<FieldA, FieldB>::PromotedType;
// We exploit that the rows of AB^T are the columns of (AB^T)^T = BA^T.
// Hence we get the row-vectors of AB^T by multiplying B to the row-vectors
// of A.
if constexpr(IsStaticSizeMatrix_v<MatrixA> and IsStaticSizeMatrix_v<WrappedMatrix>)
{
FieldMatrix<Field, MatrixA::rows, WrappedMatrix::rows> result;
for (std::size_t j=0; j<MatrixA::rows; ++j)
matrixB.wrappedMatrix().mv(matrixA[j], result[j]);
return result;
}
else
{
DynamicMatrix<Field> result(matrixA.N(), matrixB.wrappedMatrix().N());
for (std::size_t j=0; j<matrixA.N(); ++j)
matrixB.wrappedMatrix().mv(matrixA[j], result[j]);
return result;
}
}
constexpr size_type N () const
{
return wrappedMatrix().M();
}
constexpr size_type M () const
{
return wrappedMatrix().N();
}
template<class X, class Y>
void mv (const X& x, Y& y) const
{
wrappedMatrix().mtv(x,y);
}
template<class X, class Y>
void mtv (const X& x, Y& y) const
{
wrappedMatrix().mv(x,y);
}
// Return a classical representation of the matrix.
// Since we do not know the internals of the wrapped
// matrix, this will always be a dense matrix. Depending
// on whether the matrix has static size or not, this
// will be either a FieldMatrix or a DynamicMatrix.
TransposedDenseMatrixTraits_t<WrappedMatrix> asDense() const
{
TransposedDenseMatrixTraits_t<WrappedMatrix> MT;
if constexpr(not IsStaticSizeMatrix<WrappedMatrix>::value)
{
MT.resize(wrappedMatrix().M(), wrappedMatrix().N(), 0);
}
for(auto&& [M_i, i] : Dune::sparseRange(wrappedMatrix()))
for(auto&& [M_ij, j] : Dune::sparseRange(M_i))
MT[j][i] = M_ij;
return MT;
}
private:
Matrix matrix_;
};
template<class M>
using MemberFunctionTransposedConcept = std::void_t<decltype(std::declval<M>().transposed())>;
template<class M>
struct HasMemberFunctionTransposed : public Dune::Std::is_detected<MemberFunctionTransposedConcept, M> {};
} // namespace Impl
template<class Matrix,
std::enable_if_t<Impl::HasMemberFunctionTransposed<Matrix>::value, int> = 0>
auto transpose(const Matrix& matrix) {
return matrix.transposed();
}
template<class Matrix,
std::enable_if_t<not Impl::HasMemberFunctionTransposed<std::decay_t<Matrix>>::value, int> = 0>
auto transpose(Matrix&& matrix) {
return Impl::TransposedMatrixWrapper(std::forward<Matrix>(matrix));
}
template<class Matrix>
auto transpose(const std::reference_wrapper<Matrix>& matrix) {
return Impl::TransposedMatrixWrapper(matrix);
}
template<class Matrix>
auto transposedView(const Matrix& matrix) {
return transpose(std::cref(matrix));
}
} // namespace Dune
#endif // DUNE_COMMON_TRANSPOSE_HH
type_traits.hh
dynmatrix.hh
This file implements a dense matrix with dynamic numbers of rows and columns.
referencehelper.hh
boundschecking.hh
Macro for wrapping boundary checks.
matrixconcepts.hh
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
promotiontraits.hh
Compute type of the result of an arithmetic operation involving two different number types.
DUNE_ASSERT_BOUNDS
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition boundschecking.hh:30
Dune::ResolveRef_t
std::remove_reference_t< decltype(Dune::resolveRef(std::declval< T & >()))> ResolveRef_t
Type trait to resolve std::reference_wrapper.
Definition referencehelper.hh:100
Dune::sparseRange
auto sparseRange(Range &&range)
Allow structured-binding for-loops for sparse iterators.
Definition rangeutilities.hh:716
Dune::resolveRef
constexpr T & resolveRef(T &gf) noexcept
Helper function to resolve std::reference_wrapper.
Definition referencehelper.hh:47
Dune::Std::is_detected
typename detected_or< nonesuch, Op, Args... >::value_t is_detected
Detects whether Op<Args...> is valid.
Definition type_traits.hh:141
Dune
Dune namespace
Definition alignedallocator.hh:13
Dune::transpose
auto transpose(const Matrix &matrix)
Return the transposed of the given matrix.
Definition transpose.hh:234
Dune::transposedView
auto transposedView(const Matrix &matrix)
Create a view modelling the transposed matrix.
Definition transpose.hh:312
Dune::FieldTraits
Definition ftraits.hh:26
Dune::FieldTraits::field_type
T field_type
export the type representing the field
Definition ftraits.hh:28
Dune::FieldTraits::real_type
T real_type
export the type representing the real type of the field
Definition ftraits.hh:30
Dune::PromotionTraits::PromotedType
decltype(std::declval< T1 >()+std::declval< T2 >()) PromotedType
Definition promotiontraits.hh:28