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