dune-common  2.9.0
transpose.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 // SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
4 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5 #ifndef DUNE_COMMON_TRANSPOSE_HH
6 #define DUNE_COMMON_TRANSPOSE_HH
7 
8 #include <cstddef>
9 #include <functional>
10 
12 #include <dune/common/fmatrix.hh>
15 #include <dune/common/dynmatrix.hh>
17 
18 namespace Dune {
19 
20 namespace Impl {
21 
22 
23 
24  template<class M, bool = IsStaticSizeMatrix<M>::value>
25  struct TransposedDenseMatrixTraits
26  {
27  using type = Dune::FieldMatrix<typename FieldTraits<M>::field_type, M::cols, M::rows>;
28  };
29 
30  template<class M>
31  struct TransposedDenseMatrixTraits<M, false>
32  {
34  };
35 
36  template<class M>
37  using TransposedDenseMatrixTraits_t = typename TransposedDenseMatrixTraits<M>::type;
38 
39 
40 
41  // CRTP mixin class to provide the static part of the interface,
42  // namely enums rows and cols.
43  template<class WM, bool staticSize = IsStaticSizeMatrix<WM>::value>
44  class TransposedMatrixWrapperMixin {};
45 
46  template<class WM>
47  class TransposedMatrixWrapperMixin<WM, true>
48  {
49  public:
50 
52  constexpr static int rows = WM::cols;
54  constexpr static int cols = WM::rows;
55  };
56 
57 
58 
59  template<class M>
60  class TransposedMatrixWrapper;
61 
62 } // namespace Impl
63 
64 // Specialization of FieldTraits needs to be in namespace Dune::
65 template<class M>
66 struct FieldTraits< Impl::TransposedMatrixWrapper<M> >
67 {
68  using field_type = typename FieldTraits<ResolveRef_t<M>>::field_type;
69  using real_type = typename FieldTraits<ResolveRef_t<M>>::real_type;
70 };
71 
72 namespace Impl {
73 
74  // Wrapper representing the transposed of a matrix.
75  // Creating the wrapper does not compute anything
76  // but only serves for tagging the wrapped matrix
77  // for transposition. This class will store M by value.
78  // To support reference-semantic, it supports using
79  // M=std::reference_wrapper<OriginalMatrixType>.
80  template<class M>
81  class TransposedMatrixWrapper :
82  public TransposedMatrixWrapperMixin<ResolveRef_t<M>>
83  {
84  constexpr static bool hasStaticSize = IsStaticSizeMatrix<ResolveRef_t<M>>::value;
85  public:
86  using WrappedMatrix = ResolveRef_t<M>;
87 
88  const WrappedMatrix& wrappedMatrix() const {
89  return resolveRef(matrix_);
90  }
91 
92 
93  using value_type = typename WrappedMatrix::value_type;
94  using field_type = typename FieldTraits<WrappedMatrix>::field_type;
95 
96  TransposedMatrixWrapper(M&& matrix) : matrix_(std::move(matrix)) {}
97  TransposedMatrixWrapper(const M& matrix) : matrix_(matrix) {}
98 
99  template<class MatrixA,
100  std::enable_if_t<
101  ((not Impl::IsFieldMatrix<MatrixA>::value) or (not hasStaticSize))
102  and Impl::IsDenseMatrix_v<MatrixA>, int> = 0>
103  friend auto operator* (const MatrixA& matrixA, const TransposedMatrixWrapper& matrixB)
104  {
105  using FieldA = typename FieldTraits<MatrixA>::field_type;
106  using FieldB = typename FieldTraits<TransposedMatrixWrapper>::field_type;
107  using Field = typename PromotionTraits<FieldA, FieldB>::PromotedType;
108 
109  // We exploit that the rows of AB^T are the columns of (AB^T)^T = BA^T.
110  // Hence we get the row-vectors of AB^T by mutiplying B to the row-vectors
111  // of A.
112  if constexpr(IsStaticSizeMatrix_v<MatrixA> and IsStaticSizeMatrix_v<WrappedMatrix>)
113  {
114  FieldMatrix<Field, MatrixA::rows, WrappedMatrix::rows> result;
115  for (std::size_t j=0; j<MatrixA::rows; ++j)
116  matrixB.wrappedMatrix().mv(matrixA[j], result[j]);
117  return result;
118  }
119  else
120  {
121  DynamicMatrix<Field> result(matrixA.N(), matrixB.wrappedMatrix().N());
122  for (std::size_t j=0; j<matrixA.N(); ++j)
123  matrixB.wrappedMatrix().mv(matrixA[j], result[j]);
124  return result;
125  }
126  }
127 
128  template<class X, class Y>
129  void mv (const X& x, Y& y) const
130  {
131  wrappedMatrix().mtv(x,y);
132  }
133 
134  template<class X, class Y>
135  void mtv (const X& x, Y& y) const
136  {
137  wrappedMatrix().mv(x,y);
138  }
139 
140  // Return a classical representation of the matrix.
141  // Since we do not know the internals of the wrapped
142  // matrix, this will always be a dense matrix. Depending
143  // on whether the matrix has static size or not, this
144  // will be either a FieldMatrix or a DynamicMatrix.
145  TransposedDenseMatrixTraits_t<WrappedMatrix> asDense() const
146  {
147  TransposedDenseMatrixTraits_t<WrappedMatrix> MT;
148  if constexpr(not IsStaticSizeMatrix<WrappedMatrix>::value)
149  {
150  MT.resize(wrappedMatrix().M(), wrappedMatrix().N(), 0);
151  }
152  for(auto&& [M_i, i] : Dune::sparseRange(wrappedMatrix()))
153  for(auto&& [M_ij, j] : Dune::sparseRange(M_i))
154  MT[j][i] = M_ij;
155  return MT;
156  }
157 
158  private:
159 
160  M matrix_;
161  };
162 
163  template<class M>
164  using MemberFunctionTransposedConcept = std::void_t<decltype(std::declval<M>().transposed())>;
165 
166  template<class M>
167  struct HasMemberFunctionTransposed : public Dune::Std::is_detected<MemberFunctionTransposedConcept, M> {};
168 
169 } // namespace Impl
170 
171 
172 
181 template<class Matrix,
182  std::enable_if_t<Impl::HasMemberFunctionTransposed<Matrix>::value, int> = 0>
183 auto transpose(const Matrix& matrix) {
184  return matrix.transposed();
185 }
186 
187 
188 
210 template<class Matrix,
211  std::enable_if_t<not Impl::HasMemberFunctionTransposed<std::decay_t<Matrix>>::value, int> = 0>
212 auto transpose(Matrix&& matrix) {
213  return Impl::TransposedMatrixWrapper(std::forward<Matrix>(matrix));
214 }
215 
216 
217 
244 template<class Matrix>
245 auto transpose(const std::reference_wrapper<Matrix>& matrix) {
246  return Impl::TransposedMatrixWrapper(matrix);
247 }
248 
249 
250 
260 template<class Matrix>
261 auto transposedView(const Matrix& matrix) {
262  return transpose(std::cref(matrix));
263 }
264 
265 
266 
267 
268 
269 } // namespace Dune
270 
271 #endif // DUNE_COMMON_TRANSPOSE_HH
This file implements a dense matrix with dynamic numbers of rows and columns.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Compute type of the result of an arithmetic operation involving two different number types.
auto sparseRange(Range &&range)
Allow structured-binding for-loops for sparse iterators.
Definition: rangeutilities.hh:821
constexpr T & resolveRef(T &gf) noexcept
Helper function to resolve std::reference_wrapper.
Definition: referencehelper.hh:47
typename detected_or< nonesuch, Op, Args... >::value_t is_detected
Detects whether Op<Args...> is valid.
Definition: type_traits.hh:141
Dune namespace.
Definition: alignedallocator.hh:13
auto transpose(const Matrix &matrix)
Return the transposed of the given matrix.
Definition: transpose.hh:183
auto transposedView(const Matrix &matrix)
Create a view modelling the transposed matrix.
Definition: transpose.hh:261
A dense n x m matrix.
Definition: fmatrix.hh:117
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
T field_type
export the type representing the field
Definition: ftraits.hh:28
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:30
decltype(std::declval< T1 >()+std::declval< T2 >()) typedef PromotedType
Definition: promotiontraits.hh:28