Merge branch 'feature/transposed-matrix-multiply' into 'master'

Add wrapper representing the transposed of a matrix

See merge request core/dune-common!889

dune/common/CMakeLists.txt

@@ -107,6 +107,7 @@ install(FILES
         stringutility.hh
         to_unique_ptr.hh
         timer.hh
+        transpose.hh
         tupleutility.hh
         tuplevector.hh
         typelist.hh

dune/common/test/CMakeLists.txt

@@ -396,6 +396,10 @@ dune_add_test(NAME testdebugallocator_fail5
 dune_add_test(SOURCES testfloatcmp.cc
               LABELS quick)
 
+dune_add_test(SOURCES transposetest.cc
+              LINK_LIBRARIES dunecommon
+              LABELS quick)
+
 dune_add_test(SOURCES to_unique_ptrtest.cc
               LINK_LIBRARIES dunecommon
               LABELS quick)

dune/common/test/transposetest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+
+#include <type_traits>
+#include <string>
+#include <cstddef>
+
+#include <dune/common/transpose.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/classname.hh>
+#include <dune/common/test/testsuite.hh>
+
+
+
+// check A*transpose(B)
+template<class A, class B>
+auto checkAxBT(const A& a, const B&b)
+{
+  Dune::TestSuite suite(std::string{"Check transpose with A="} + Dune::className<A>() + " and B=" + Dune::className<B>());
+
+  // compute abt
+  auto abt = a * transpose(b);
+
+  // check result type
+  using FieldA = typename Dune::FieldTraits<A>::field_type;
+  using FieldB = typename Dune::FieldTraits<B>::field_type;
+  using Field = typename Dune::PromotionTraits<FieldA, FieldB>::PromotedType;
+  using ABT = Dune::FieldMatrix<Field,A::rows, B::rows>;
+  suite.check(std::is_same<decltype(abt), ABT>())
+    << "Result type of A*transpose(B) should be " << Dune::className<ABT>() << " but is " << Dune::className<decltype(abt)>();
+
+  // manually compute result value
+  auto abt_check = ABT{};
+  for(std::size_t i=0; i<A::rows; ++i)
+    for(std::size_t j=0; j<B::rows; ++j)
+      for(auto&& [b_jk, k] : Dune::sparseRange(b[j]))
+        abt_check[i][j] += a[i][k]*b_jk;
+
+  // check result value
+  bool equal = true;
+  for(std::size_t i=0; i<A::rows; ++i)
+    for(std::size_t j=0; j<B::rows; ++j)
+      equal = equal and (abt_check[i][j] == abt[i][j]);
+  suite.check(equal)
+    << "Result of A*transpose(B) should be \n" << abt_check << " but is \n" << abt;
+
+  return suite;
+}
+
+
+int main()
+{
+  Dune::TestSuite suite;
+
+  // fill dense matrix with test data
+  auto testFillDense = [](auto& matrix) {
+    std::size_t k=0;
+    for(std::size_t i=0; i<matrix.N(); ++i)
+      for(std::size_t j=0; j<matrix.M(); ++j)
+        matrix[i][j] = k++;
+  };
+
+  {
+    auto a = Dune::FieldMatrix<double,3,4>{};
+    auto b = Dune::FieldMatrix<double,7,4>{};
+    testFillDense(a);
+    testFillDense(b);
+    suite.subTest(checkAxBT(a,b));
+  }
+
+  {
+    auto a = Dune::FieldMatrix<double,1,2>{};
+    auto b = Dune::FieldMatrix<double,3,2>{};
+    testFillDense(a);
+    testFillDense(b);
+    suite.subTest(checkAxBT(a,b));
+  }
+
+  {
+    auto a = Dune::FieldMatrix<double,1,2>{};
+    auto b = Dune::FieldMatrix<double,1,2>{};
+    testFillDense(a);
+    testFillDense(b);
+    suite.subTest(checkAxBT(a,b));
+  }
+
+  {
+    auto a = Dune::FieldMatrix<double,3,4>{};
+    auto b = Dune::DiagonalMatrix<double,4>{};
+    testFillDense(a);
+    b = {0, 1, 2, 3};
+    suite.subTest(checkAxBT(a,b));
+  }
+
+  return suite.exit();
+}

dune/common/transpose.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_COMMON_TRANSPOSE_HH
+#define DUNE_COMMON_TRANSPOSE_HH
+
+#include <cstddef>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/promotiontraits.hh>
+
+namespace Dune {
+
+namespace Impl {
+
+  // Wrapper representing the transposed of a matrix.
+  // Creating the wrapper does not compute anything
+  // but only serves for tagging the wrapped matrix
+  // for transposition.
+  template<class M>
+  class TransposedMatrixWrapper
+  {
+  public:
+
+    enum {
+      //! The number of rows.
+      rows = M::cols,
+      //! The number of columns.
+      cols = M::rows
+    };
+
+    TransposedMatrixWrapper(const M& matrix) : matrix_(matrix) {}
+    TransposedMatrixWrapper(const TransposedMatrixWrapper&) = delete;
+    TransposedMatrixWrapper(TransposedMatrixWrapper&&) = delete;
+
+    template<class OtherField, int otherRows>
+    friend auto operator* (const FieldMatrix<OtherField, otherRows, rows>& matrixA,
+                            const TransposedMatrixWrapper& matrixB)
+    {
+      using ThisField = typename FieldTraits<M>::field_type;
+      using Field = typename PromotionTraits<ThisField, OtherField>::PromotedType;
+      FieldMatrix<Field, otherRows, cols> result;
+      for (std::size_t j=0; j<otherRows; ++j)
+        matrixB.matrix_.mv(matrixA[j], result[j]);
+      return result;
+    }
+
+  private:
+
+    const M& matrix_;
+  };
+
+} // namespace Impl
+
+/**
+ * \brief Create a wrapper modelling the transposed matrix
+ *
+ * Currently the wrapper only implements
+ * \code
+ * auto c = a*transpose(b);
+ * \endcode
+ * if a is a FieldMatrix of appropriate size. This is
+ * optimal even for sparse b because it only relies on
+ * calling b.mv(a[i], c[i]) for the rows of a.
+ *
+ * Since the created object only stores a reference
+ * to the wrapped matrix, it cannot be modified and
+ * should not be stored but used directly.
+ */
+template<class Matrix>
+auto transpose(const Matrix& matrix) {
+  return Impl::TransposedMatrixWrapper<Matrix>(matrix);
+}
+
+
+} // namespace Dune
+
+#endif // DUNE_COMMON_TRANSPOSE_HH