Merge branch 'feature/cleanup-istlmatrixbackend' into 'master'

[cleanup]Simplify ISTLMatrixBackend entry access by multi-indices

See merge request !149

dune/fufem/backends/istlmatrixbackend.hh

@@ -23,9 +23,9 @@ namespace Dune::Fufem {
 
 
 
-template<class M, class E=typename M::field_type>
-class ISTLMatrixBackend
-{
+namespace Impl {
+
+
 
   // The following is essentially a copy of Dune::Functions::hybridIndexAccess
   // But since matrices do not provide size(), we have to use N() instead.
@@ -49,92 +49,112 @@ class ISTLMatrixBackend
         });
   }
 
-  template<class Result, class RowIndex, class ColIndex>
-  struct MatrixMultiIndexResolver
+
+
+  // Call f(matrix[i][j]) and return the result.
+  // This works with dynamic indices i,j, even for a multi-type matrix.
+  // However, it requires that f() has a unique return type.
+  template<class Matrix, class RowIndex, class ColIndex, class F>
+  decltype(auto) visitMatrixEntry(Matrix&& matrix, const RowIndex& i, const ColIndex& j, F&& f)
   {
-    MatrixMultiIndexResolver(const RowIndex& rowIndex, const ColIndex& colIndex) :
-      rowIndex_(rowIndex),
-      colIndex_(colIndex)
-    {}
-
-    template<class C>
-    using isReturnable =
-    std::integral_constant<bool,
-      std::is_convertible<C&, Result>::value or
-      Dune::MatrixVector::Traits::ScalarTraits<std::decay_t<C>>::isScalar
-    >;
-
-    template<class C>
-    auto& toScalar(C&& c) {
-      return std::forward<C>(c);
-    }
+    using namespace Dune::Functions;
+    return hybridRowIndexAccess(matrix, i, [&](auto&& matrix_i) -> decltype(auto) {
+      return hybridIndexAccess(matrix_i, j, f);
+    });
+  }
 
-    // Catch the cases where a FieldMatrix<K, 1, 1> was supplied
-    auto& toScalar(FieldMatrix<std::decay_t<Result>, 1, 1>& c) {
-      return c[0][0];
-    }
 
-    auto& toScalar(const FieldMatrix<std::decay_t<Result>, 1, 1>& c) {
-      return c[0][0];
-    }
 
-    template<class C,
-      typename std::enable_if<not isReturnable<C>::value, int>::type = 0,
-      typename std::enable_if<IsSingleRowMatrix<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+  // Call f(matrix[i0][j0]...[in][jm]), by recursively resolving row- and column-multi-indices.
+  // Whenenver a SingleRowMatrix or SingleColumnMatrix is encountered, a zero row- or column-index
+  // is inserted, respectively. The recursion ends, whenever f(matrixEntry) can be called.
+  // This works with dynamic indices i,j, even for a multi-type matrix.
+  // However, it requires that f() has a unique return type.
+  template<class Matrix, class RowIndex, class ColIndex, class F>
+  static decltype(auto) visitMatrixEntryRecursive(Matrix& matrix, const RowIndex& iii, const ColIndex& jjj, F&& f)
+  {
+    using namespace Dune::Indices;
+    using namespace Dune::Functions::Imp;
+    static constexpr bool isSingleRow = IsSingleRowMatrix<Matrix>::value;
+    static constexpr bool isSingleCol = IsSingleColumnMatrix<Matrix>::value;
+    auto splitIndex = [] (auto&& multiIndex) { return std::make_tuple(multiIndex[_0], shiftedDynamicMultiIndex<1>(multiIndex)) ;};
+    if constexpr (std::is_invocable_v<F, Matrix&>)
     {
-      using namespace Dune::Indices;
-      using namespace Dune::Functions;
-      auto&& subRowIndex = rowIndex_;
-      auto&& subColIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(colIndex_);
-      auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return static_cast<Result>((hybridRowIndexAccess(c, _0, [&](auto&& ci) -> decltype(auto) {
-        return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
-      })));
+      return f(matrix);
     }
-
-    template<class C,
-      typename std::enable_if<not isReturnable<C>::value, int>::type = 0,
-      typename std::enable_if<IsSingleColumnMatrix<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+    else if constexpr (isSingleRow)
     {
-      using namespace Dune::Indices;
-      using namespace Dune::Functions;
-      auto&& subRowIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(rowIndex_);
-      auto&& subColIndex = colIndex_;
-      auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return static_cast<Result>((hybridRowIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
-        return static_cast<Result>((hybridIndexAccess(ci, _0, subIndexResolver)));
-      })));
+      assert(jjj.size()>0);
+      auto [j, jj] = splitIndex(jjj);
+      // We have to capture jj explicitly by value, because capturing structured bindings
+      // by reference is not allowed before C++20
+      return visitMatrixEntry(matrix, _0, j, [&, jj=jj](auto&& matrix_0j) -> decltype(auto) {
+        return visitMatrixEntryRecursive(matrix_0j, iii, jj, f);
+      });
     }
-
-    template<class C,
-      typename std::enable_if<not isReturnable<C>::value, int>::type = 0,
-      typename std::enable_if<not IsSingleRowMatrix<C>::value, int>::type = 0,
-      typename std::enable_if<not IsSingleColumnMatrix<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+    else if constexpr (isSingleCol)
     {
-      using namespace Dune::Indices;
-      using namespace Dune::Functions;
-      auto&& subRowIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(rowIndex_);
-      auto&& subColIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(colIndex_);
-      auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
-      return static_cast<Result>((hybridRowIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
-        return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
-      })));
+      assert(iii.size()>0);
+      auto [i, ii] = splitIndex(iii);
+      // We have to capture ii explicitly by value, because capturing structured bindings
+      // by reference is not allowed before C++20
+      return visitMatrixEntry(matrix, i, _0, [&, ii=ii](auto&& matrix_i0) -> decltype(auto) {
+        return visitMatrixEntryRecursive(matrix_i0, ii, jjj, f);
+      });
     }
-
-    template<class C,
-      typename std::enable_if<isReturnable<C>::value, int>::type = 0>
-    Result operator()(C&& c)
+    else
     {
-      return static_cast<Result>(toScalar(std::forward<C>(c)));
+      assert(iii.size()>0);
+      assert(jjj.size()>0);
+      auto [i, ii] = splitIndex(iii);
+      auto [j, jj] = splitIndex(jjj);
+      // We have to capture ii and jj explicitly by value, because capturing structured bindings
+      // by reference is not allowed before C++20
+      return visitMatrixEntry(matrix, i, j, [&, ii=ii, jj=jj](auto&& matrix_ij) -> decltype(auto) {
+        return visitMatrixEntryRecursive(matrix_ij, ii, jj, f);
+      });
+    }
+  }
+
+
+
+} // namespace Impl
+
+
+
+template<class M, class E=typename M::field_type>
+class ISTLMatrixBackend
+{
+
+  template<class Result>
+  struct ToScalar;
+
+  template<class R>
+  struct ToScalar<R&>
+  {
+    template<class Matrix,
+      std::enable_if_t<std::is_convertible_v<Matrix&, R&>, int> = 0>
+    R& operator()(Matrix& matrix) {
+      return matrix;
     }
 
+    R& operator()(Dune::FieldMatrix<R, 1, 1>& matrix) {
+      return matrix[0][0];
+    }
+  };
 
+  template<class R>
+  struct ToScalar<const R&>
+  {
+    template<class Matrix,
+      std::enable_if_t<std::is_convertible_v<Matrix&, const R&>, int> = 0>
+    const R& operator()(Matrix& matrix) {
+      return matrix;
+    }
 
-    const RowIndex& rowIndex_;
-    const ColIndex& colIndex_;
+    const R& operator()(const Dune::FieldMatrix<R, 1, 1>& matrix) {
+      return matrix[0][0];
+    }
   };
 
 
@@ -157,15 +177,13 @@ public:
   template<class RowMultiIndex, class ColMultiIndex>
   const Entry& operator()(const RowMultiIndex& row, const ColMultiIndex& col) const
   {
-    MatrixMultiIndexResolver<const Entry&, RowMultiIndex, ColMultiIndex> i(row, col);
-    return i(*matrix_);
+    return Impl::visitMatrixEntryRecursive(*matrix_, row, col, ToScalar<const Entry&>());
   }
 
   template<class RowMultiIndex, class ColMultiIndex>
   Entry& operator()(const RowMultiIndex& row, const ColMultiIndex& col)
   {
-    MatrixMultiIndexResolver<Entry&, RowMultiIndex, ColMultiIndex> i(row, col);
-    return i(*matrix_);
+    return Impl::visitMatrixEntryRecursive(*matrix_, row, col, ToScalar<Entry&>());
   }
 
   /**

dune/fufem/test/istlbackendtest.cc

 #include <config.h>
+
 #include <dune/common/test/testsuite.hh>
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/indices.hh>
+#include <dune/common/rangeutilities.hh>
+#include <dune/common/classname.hh>
 
 #include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/multitypeblockmatrix.hh>
 #include <dune/istl/multitypeblockvector.hh>
 
-#include <dune/fufem/backends/istlmatrixbackend.hh>
-#include <dune/fufem/test/common.hh>
-
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/functionspacebases/powerbasis.hh>
 #include <dune/functions/functionspacebases/compositebasis.hh>
-#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/backends/istlmatrixbackend.hh>
+#include <dune/fufem/test/common.hh>
 
 // Due to a formerly unnoticed bug, activate bounds checking:
 #define DUNE_CHECK_BOUNDS 1
 
-using namespace Dune;
 
-template<class Entry, class Matrix, class RowIndex, class ColIndex>
-TestSuite testISTLBackend(Matrix& mat, RowIndex&& row, ColIndex&& col) {
-  TestSuite suite;
 
-  auto backend = Fufem::istlMatrixBackend<Entry>(mat);
+// Build matrix that contains entries for all pairs (i,j)
+// of basis functions that share an element.
+template<class Matrix, class Basis>
+Matrix buildMatrixForBasis(const Basis& basis) {
+  Matrix matrix;
 
-  suite.check(backend(row, col) == 1.0, "Check matrix access") << "Type" << typeid(mat).name() << " failed";
+  auto backend = Dune::Fufem::istlMatrixBackend(matrix);
+  auto pattern = backend.patternBuilder();
 
-  return suite;
+  pattern.resize(basis,basis);
+  auto localView = basis.localView();
+  for(const auto& e : elements(basis.gridView()))
+  {
+    localView.bind(e);
+    for(auto i : Dune::range(localView.size()))
+      for(auto j : Dune::range(localView.size()))
+        pattern.insertEntry(localView.index(i), localView.index(j));
+  }
+  pattern.setupMatrix();
+  return matrix;
 }
 
-template<class Matrix>
-TestSuite testMultiType(Matrix& mat) {
-  TestSuite suite;
 
-  auto backend = Fufem::istlMatrixBackend(mat);
-  auto pattern = backend.patternBuilder();
 
-  // create small grid
-  auto grid = constructCoarseYaspGrid<3>();
-  auto gridView = grid->leafGridView();
+// Test istlMatrixBackend() for existing matrix.
+template<class Entry, class Matrix, class Basis>
+Dune::TestSuite testISTLMatrixBackendForBasis(Matrix& matrix, const Basis& basis) {
+  Dune::TestSuite suite(Dune::className(matrix));
 
-  // create small basis
-  namespace BB = Functions::BasisBuilder;
-
-  // create a composite basis matching the matrix type
-  const auto basis = BB::makeBasis(
-                        gridView,
-                        BB::composite(
-                          BB::power<3>(
-                            BB::lagrange<1>()
-                          ),
-                          BB::power<5>(
-                            BB::lagrange<0>()
-                          )
-                        )
-                      );
+// FieldTraits<MultiTypeBlockMatrix> are not implemented properly
+// in dune-2.9. Hence determining the field type automatically
+// in istlMatrixBackend fails and we have to provide the type
+// manually.
+#if DUNE_VERSION_GT(DUNE_COMMON, 2, 9)
+  auto backend = Dune::Fufem::istlMatrixBackend(matrix);
+#else
+  auto backend = Dune::Fufem::istlMatrixBackend<Entry>(matrix);
+#endif
 
-  pattern.resize(basis,basis);
+  backend.assign(1.0);
 
-  // set the active entries
   auto localView = basis.localView();
-  for(const auto& e : elements(gridView))
+  for(const auto& e : elements(basis.gridView()))
   {
     localView.bind(e);
-    for(size_t i=0; i<localView.size(); i++)
+    for(auto i : Dune::range(localView.size()))
     {
-      auto idx = localView.index(i);
-      pattern.insertEntry(idx,idx);
+      auto rowIndex = localView.index(i);
+      for(auto j : Dune::range(localView.size()))
+      {
+        auto colIndex = localView.index(j);
+        suite.check(backend(rowIndex, colIndex) == 1.0, "Multi-index access")
+          << "Unexpected entry (" << rowIndex << "," << colIndex << ")";
+      }
     }
   }
 
-  pattern.setupMatrix();
-  mat = 1.0;
+  return suite;
+}
 
-  for(const auto& e : elements(gridView))
-  {
-    localView.bind(e);
-    for(size_t i=0; i<localView.size(); i++)
-    {
-      auto idx = localView.index(i);
-      suite.check(backend(idx, idx) == 1.0,
-                  "Check matrix access") << "Type" << typeid(mat).name() << " failed";
-    }
-  }
 
-  return suite;
+
+// Build matrix for a basis and test istlMatrixBackend() for it.
+template<class Entry, class Matrix, class Basis>
+Dune::TestSuite testISTLMatrixBackendForBasis(const Basis& basis) {
+  Matrix matrix = buildMatrixForBasis<Matrix>(basis);
+  return testISTLMatrixBackendForBasis<Entry>(matrix, basis);
 }
 
-int main(int argc, char** argv) {
-  MPIHelper::instance(argc, argv);
 
-  TestSuite suite;
 
-  using K = double;
+template<class... V>
+using MTBV = Dune::MultiTypeBlockVector<V...>;
 
-  // Matrix<FieldMatrix<K, 1,1>
-  {
-    using Matrix = Matrix<FieldMatrix<K,1,1>>;
-    Matrix mat(2,2);
-    mat =1.0;
+template<class... V>
+using MTBM = Dune::MultiTypeBlockMatrix<V...>;
 
-    // check with multiindex with length 1
-    {
-      auto idx = std::array<size_t, 1>{{1}};
+template<class K, int n, int m>
+using FM = Dune::FieldMatrix<K,n,m>;
 
-      suite.subTest(testISTLBackend<K>(mat, idx, idx));
+template<class M>
+using BCRSM = Dune::BCRSMatrix<M>;
 
-      // check const:
-      const auto& const_mat = mat;
-      suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
-    }
-    // check with multiindex with length 2
-    {
-      // someone might set a trailing zero to account for the matrix
-      // character of FieldMatrix<K, 1,1>:
-      auto idx = std::array<size_t, 2>{{1, 0}};
+template<class M>
+using SRM = Dune::Fufem::SingleRowMatrix<M>;
 
-      suite.subTest(testISTLBackend<K>(mat, idx, idx));
+template<class M>
+using SCM = Dune::Fufem::SingleColumnMatrix<M>;
 
-      // check const:
-      const auto& const_mat = mat;
-      suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
-    }
+
+
+int main(int argc, char** argv) {
+  Dune::MPIHelper::instance(argc, argv);
+
+  Dune::TestSuite suite;
+
+  using K = double;
+  using namespace Dune::Functions::BasisBuilder;
+
+  // Construct grid for checking backend against bases
+  auto grid = constructCoarseYaspGrid<3>();
+  auto gridView = grid->leafGridView();
+
+  //****************************************************************************
+  // Test with dense matrix and scalar entries
+  //****************************************************************************
+  {
+    using Matrix = Dune::Matrix<K>;
+    const auto basis = makeBasis(
+                          gridView,
+                          lagrange<1>()
+                        );
+    // MatrixBuilder is not specialized for this type yet.
+    Matrix matrix(basis.size(), basis.size());
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
   }
 
-  // Matrix<FieldMatrix<K, n,n> n>1
+  //****************************************************************************
+  // Test with dense matrix and scalar entries (as FM<K,1,1>
+  //****************************************************************************
   {
-    constexpr int n = 2;
-    using Matrix = Matrix<FieldMatrix<K, n, n>>;
-    Matrix mat(2,2);
-    mat = 1.0;
+    using Matrix = Dune::Matrix<FM<K,1,1>>;
+    const auto basis = makeBasis(
+                          gridView,
+                          lagrange<1>()
+                        );
+    Matrix matrix(basis.size(), basis.size());
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
+  }
+
+  //****************************************************************************
+  // Test with sparse matrix and scalar entries
+  //****************************************************************************
+  {
+    using Matrix = BCRSM<K>;
+    const auto basis = makeBasis(
+                          gridView,
+                          lagrange<1>()
+                        );
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
+  }
 
-    auto idx = std::array<size_t, 2>{{1, 1}};
+  //****************************************************************************
+  // Test with sparse matrix and scalar entries (as FM<K,1,1>)
+  //****************************************************************************
+  {
+    using Matrix = BCRSM<FM<K,1,1>>;
+    const auto basis = makeBasis(
+                          gridView,
+                          lagrange<1>()
+                        );
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
+  }
 
-    suite.subTest(testISTLBackend<K>(mat, idx, idx));
+  //****************************************************************************
+  // Test with dense matrix and block entries
+  //****************************************************************************
+  {
+    using Matrix = Dune::Matrix<FM<K,2,2>>;
+    const auto basis = makeBasis(
+                          gridView,
+                          power<2>(
+                            lagrange<1>()
+                          )
+                        );
+    // MatrixBuilder is not specialized for this type yet.
+    Matrix matrix(basis.size(), basis.size());
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
+  }
 
-    // check const:
-    const auto& const_mat = mat;
-    suite.subTest(testISTLBackend<const K>(const_mat, idx, idx));
+  //****************************************************************************
+  // Test with dense matrix and dynamic block entries
+  //****************************************************************************
+  {
+    using Matrix = Dune::Matrix<Dune::Matrix<K>>;
+    const auto basis = makeBasis(
+                          gridView,
+                          power<2>(
+                            lagrange<1>()
+                          )
+                        );
+    // MatrixBuilder is not specialized for this type yet.
+    Matrix matrix(basis.size(), basis.size());
+    for (auto i: Dune::range(basis.size()))
+      for (auto j: Dune::range(basis.size()))
+        matrix[i][j].setSize(basis.size({i}), basis.size({j}));
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
+  }
+  //****************************************************************************
+  // Test with sparse matrix and block entries
+  //****************************************************************************
+  {
+    using Matrix = BCRSM<FM<K,2,2>>;
+    const auto basis = makeBasis(
+                          gridView,
+                          power<2>(
+                            lagrange<1>()
+                          )
+                        );
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
   }
 
-  // MultiTypeMatrix
+  //****************************************************************************
+  // Test with MTBM and uniform indices
+  //****************************************************************************
   {
-    using Matrix = MultiTypeBlockMatrix<MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,3,3>>,BCRSMatrix<FieldMatrix<double,3,5>>>,
-                                        MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,5,3>>,BCRSMatrix<FieldMatrix<double,5,5>>>>;
+    using Matrix =
+      MTBM<
+        MTBV<  BCRSM<FM<K,3,3>>, BCRSM<FM<K,3,5>>  >,
+        MTBV<  BCRSM<FM<K,5,3>>, BCRSM<FM<K,5,5>>  >
+      >;
+    const auto basis = makeBasis(
+                          gridView,
+                          composite(
+                            power<3>(
+                              lagrange<1>()
+                            ),
+                            power<5>(
+                              lagrange<0>()
+                            )
+                          )
+                        );
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
+  }
 
-    Matrix mat;
+  //****************************************************************************
+  // Test with MTBM and non-uniform indices
+  //****************************************************************************
+  {
+    using Matrix =
+      MTBM<
+        MTBV<  BCRSM<FM<K,3,3>>,      BCRSM<SCM<FM<K,3,1>>>  >,
+        MTBV<  BCRSM<SRM<FM<K,3,1>>>, BCRSM<FM<K,1,1>>       >
+      >;
+    const auto basis = makeBasis(
+                          gridView,
+                          composite(
+                            power<3>(
+                              lagrange<2>()
+                            ),
+                            lagrange<1>()
+                          )
+                        );
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(basis));
+  }
 
-    suite.subTest(testMultiType(mat));
+  //****************************************************************************
+  // Test with nested MTBM and non-uniform indices (scalar entries as FM<K,1,1>)
+  //****************************************************************************
+  {
+    using Matrix =
+      MTBM<
+        MTBV<
+          MTBM<
+            MTBV<  Dune::Matrix<FM<K,3,3>>,      Dune::Matrix<SCM<FM<K,3,1>>>  >,
+            MTBV<  Dune::Matrix<SRM<FM<K,3,1>>>, Dune::Matrix<FM<K,1,1>>       >
+          >
+        >
+      >;
+    const auto basis = makeBasis(
+                          gridView,
+                          composite(
+                            composite(
+                              power<3>(
+                                lagrange<2>()
+                              ),
+                              lagrange<1>()
+                            )
+                          )
+                        );
+    // MatrixBuilder is not specialized for this type yet.
+    using namespace Dune::Indices;
+    Matrix matrix;
+    matrix[_0][_0][_0][_0].setSize(basis.size({0,0}), basis.size({0,0}));
+    matrix[_0][_0][_0][_1].setSize(basis.size({0,0}), basis.size({0,1}));
+    matrix[_0][_0][_1][_0].setSize(basis.size({0,1}), basis.size({0,0}));
+    matrix[_0][_0][_1][_1].setSize(basis.size({0,1}), basis.size({0,1}));
+    suite.subTest(testISTLMatrixBackendForBasis<K, Matrix>(matrix, basis));
   }
 
   return suite.exit();