DUNE PDELab (git)

transpose.hh
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 © 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
11#include <dune/common/std/type_traits.hh>
14#include <dune/common/referencehelper.hh>
16#include <dune/common/matrixconcepts.hh>
17
18namespace Dune {
19
20namespace 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::
65template<class M>
66struct 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
72namespace 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 multiplying 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
181template<class Matrix,
182 std::enable_if_t<Impl::HasMemberFunctionTransposed<Matrix>::value, int> = 0>
183auto transpose(const Matrix& matrix) {
184 return matrix.transposed();
185}
186
187
188
210template<class Matrix,
211 std::enable_if_t<not Impl::HasMemberFunctionTransposed<std::decay_t<Matrix>>::value, int> = 0>
212auto transpose(Matrix&& matrix) {
213 return Impl::TransposedMatrixWrapper(std::forward<Matrix>(matrix));
214}
215
216
217
244template<class Matrix>
245auto transpose(const std::reference_wrapper<Matrix>& matrix) {
246 return Impl::TransposedMatrixWrapper(matrix);
247}
248
249
250
260template<class Matrix>
261auto 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
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
A dense n x m matrix.
Definition: fmatrix.hh:117
A generic dynamic dense matrix.
Definition: matrix.hh:561
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 ...
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:145
auto sparseRange(Range &&range)
Allow structured-binding for-loops for sparse iterators.
Definition: rangeutilities.hh:722
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
STL namespace.
Compute type of the result of an arithmetic operation involving two different number types.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)