transpose.hh

// -*- 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