istlmatrixbackend.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_ASSEMBLER_BACKENDS_ISTLMATRIXBACKEND_HH
#define DUNE_ASSEMBLER_BACKENDS_ISTLMATRIXBACKEND_HH
 
#include <array>
#include <vector>
 
#include <dune/functions/common/indexaccess.hh>
#include <dune/functions/common/multiindex.hh>
 
#include <dune/common/rangeutilities.hh>
 
#include <dune/assembler/utility.hh>
#include <dune/assembler/backends/flatmatrix.hh>
#include <dune/assembler/backends/localpattern.hh>
#include <dune/assembler/backends/singlerowcolmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
 
 
 
namespace Dune::Assembler::Impl {
 
 
 
// Mix-in class for MatrixPatternBackend implementations.
// This CRTP base class implements scatter() by adding individual
// index pairs to the global pattern using the
// add(rowIndex, colIndex) provided by the derived
// implementation class.
// The used LocalPattern class can be switched using a
// template parameter, using  LocalPatternMatrix as default.
template <class Impl, class LP=Dune::Assembler::LocalPatternMatrix>
class MatrixPatternBackendMixIn
{
public:
  using LocalPattern = LP;
 
  // use a given pattern
  template <class RowLocalView, class ColLocalView>
  void scatter (const RowLocalView& rowLocalView, const ColLocalView& colLocalView,
                const LocalPattern& localPattern)
  {
    auto rowIndices = Dune::Assembler::Impl::localIndices(rowLocalView.tree());
    auto colIndices = Dune::Assembler::Impl::localIndices(colLocalView.tree());
    localPattern.visitEntries(rowIndices, colIndices, [&](const auto& i, const auto&j) {
      static_cast<Impl*>(this)->add(rowLocalView.index(i), colLocalView.index(j));
    });
  }
 
};
 
 
 
// Number of matrix rows provided as Dune::index_constant
// if matrix does not support dynamic access.
template <class Matrix>
auto numRows (const Matrix& m)
{
  if constexpr (requires { m[0]; })
    return m.N();
  else
    return Dune::index_constant<Matrix::N()>();
}
 
// Number of matrix columns provided as Dune::index_constant
// if matrix does not support dynamic access.
template <class Matrix>
auto numCols (const Matrix& m)
{
  if constexpr (requires { m[Dune::Indices::_0][0]; })
    return m.M();
  else
    return Dune::index_constant<Matrix::M()>();
}
 
// Invoke callback with index corresponding to i in given range
template <class T, T to, T from, class I, class F>
static decltype(auto) indexInRange(const Dune::StaticIntegralRange<T, to, from>& range, const I& i, F&& f)
{
  if constexpr (to==from+1)
    return f(std::integral_constant<T, from>());
  else
  {
    if (i==from)
      return f(std::integral_constant<T, from>());
    else
      return indexInRange(Dune::StaticIntegralRange<T, to, from+1>(), i, f);
  }
}
 
// Invoke callback with index corresponding to i in given range
template <class T, class I, class F>
static decltype(auto) indexInRange(const Dune::IntegralRange<T>& range, const I& i, F&& f)
{
  return f(i);
}
 
// This calls f(matrix[i][j]) and returns 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) hybridFlatMatrixEntryAccess(Matrix&& matrix, const RowIndex& i, const ColIndex& j, F&& f)
{
  return indexInRange(range(numRows(matrix)), i, [&](auto ii) -> decltype(auto) {
      return indexInRange(range(numCols(matrix)), j, [&](auto jj) -> decltype(auto) {
        return f(matrix[ii][jj]);
      });
    });
}
 
 
 
// 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) hybridRecursiveMatrixEntryAccess(Matrix& matrix, const RowIndex& iii, const ColIndex& jjj, F&& f)
{
  using Dune::Functions::Imp::shiftedDynamicMultiIndex;
  auto splitIndex = [] (auto&& multiIndex, auto insertZero) {
    if constexpr (decltype(insertZero)::value)
      return std::make_tuple(Dune::Indices::_0, shiftedDynamicMultiIndex<0>(multiIndex));
    else
    {
      assert(multiIndex.size()>0);
      return std::make_tuple(multiIndex[Dune::Indices::_0], shiftedDynamicMultiIndex<1>(multiIndex)) ;
    }
  };
  if constexpr (std::is_invocable_v<F, Matrix&>)
    return f(matrix);
  else if constexpr (std::is_invocable_v<F, Matrix&, const RowIndex&, const ColIndex&>)
    return f(matrix, iii, jjj);
  else
  {
    auto [i, ii] = splitIndex(iii, IsSingleRowMatrix<Matrix>());
    auto [j, jj] = splitIndex(jjj, IsSingleColumnMatrix<Matrix>());
    return hybridFlatMatrixEntryAccess(matrix, i, j, [&,ii=ii,jj=jj](auto&& matrix_ij) -> decltype(auto) {
        return hybridRecursiveMatrixEntryAccess(matrix_ij, ii, jj, f);
      });
  }
}
 
 
 
} // end namespace Dune::Assembler::Impl
 
 
 
 
namespace Dune::Assembler {
 
// forward declaration
template <class Matrix>
class ISTLMatrixPatternBackend;
 
template <class Matrix>
class ISTLMatrixBackend
{
 
public:
  using value_type = typename Matrix::field_type;
  using LocalMatrix = FlatMatrix<value_type>;
 
  using PatternBuilder = ISTLMatrixPatternBackend<Matrix>;
  using LocalPattern = typename PatternBuilder::LocalPattern;
 
private:
 
  template <class R>
  struct ConvertToRef
  {
 
    template <class M>
    R& operator() (M& matrix) const
    requires(std::is_convertible_v<M&, R&>)
    {
      return matrix;
    }
 
    // This is needed for treating `Dune::FieldMatrix<R,1,1>` like
    // a scalar, because there's no cast to `R&`. We don't need
    // an overload for the const case, because there is a cast
    // to `const R&`.
    R& operator() (Dune::FieldMatrix<R, 1, 1>& matrix) const
    requires(not std::is_const_v<R>)
    {
      return matrix[0][0];
    }
  };
 
  class AccumulateValue
  {
    const value_type& value_;
 
  public:
 
    AccumulateValue (const value_type& value)
      : value_(value)
    {}
 
    template <class M>
    void operator() (M& matrix) const
    requires(std::is_convertible_v<M&, value_type&>)
    {
      matrix += value_;
    }
 
    template <class T, class A, class RowIndex, class ColIndex>
    void operator() (Dune::BCRSMatrix<T, A>& matrix, const RowIndex& ii, const ColIndex& jj) const
    requires(std::is_convertible_v<T&, value_type&>)
    {
      using namespace Dune::Indices;
      matrix[ii[_0]][jj[_0]] += value_;
    }
 
    template <class T, int n, int m, class A, class RowIndex, class ColIndex>
    void operator() (Dune::BCRSMatrix<Dune::FieldMatrix<T, n, m>, A>& matrix, const RowIndex& ii, const ColIndex& jj) const
    requires(std::is_convertible_v<T&, value_type&>)
    {
      using namespace Dune::Indices;
      matrix[ii[_0]][jj[_0]][ii[_1]][jj[_1]] += value_;
    }
 
    template <class T, int n, class A, class RowIndex, class ColIndex>
    void operator() (Dune::BCRSMatrix<Dune::DiagonalMatrix<T, n>, A>& matrix, const RowIndex& ii, const ColIndex& jj) const
    requires(std::is_convertible_v<T&, value_type&>)
    {
      using namespace Dune::Indices;
      if (ii[_1] == jj[_1])
        matrix[ii[_0]][jj[_0]].diagonal(ii[_1]) += value_;
    }
 
    template <class T, int n, class A, class RowIndex, class ColIndex>
    void operator() (Dune::BCRSMatrix<Dune::ScaledIdentityMatrix<T, n>, A>& matrix, const RowIndex& ii, const ColIndex& jj) const
    requires(std::is_convertible_v<T&, value_type&>)
    {
      using namespace Dune::Indices;
      if ((ii[_1] == 0) and (jj[_1] == 0))
        matrix[ii[_0]][jj[_0]].scalar() += value_;
    }
 
    template <class F, int n, class RowIndex, class ColIndex>
    void operator() (Dune::DiagonalMatrix<F, n>& matrix, const RowIndex& ii, const ColIndex& jj) const
    {
      using namespace Dune::Indices;
      if (ii[_0] == jj[_0])
        matrix.diagonal(ii[_0]) += value_;
    }
 
    template <class F, int n, class RowIndex, class ColIndex>
    void operator() (Dune::ScaledIdentityMatrix<F, n>& matrix, const RowIndex& ii, const ColIndex& jj) const
    {
      using namespace Dune::Indices;
      if ((ii[_0] == 0) and (jj[_0] == 0))
        matrix.scalar() += value_;
    }
 
  };
 
public:
  ISTLMatrixBackend (Matrix& matrix)
    : matrix_(matrix)
  {}
 
  PatternBuilder patternBuilder ()
  {
    return PatternBuilder{matrix_};
  }
 
  template <class RowMultiIndex, class ColMultiIndex>
  const value_type& operator() (const RowMultiIndex& row, const ColMultiIndex& col) const
  {
    return Impl::hybridRecursiveMatrixEntryAccess(matrix_, row, col, ConvertToRef<const value_type>());
  }
 
  template <class RowMultiIndex, class ColMultiIndex>
  value_type& operator() (const RowMultiIndex& row, const ColMultiIndex& col)
  {
    return Impl::hybridRecursiveMatrixEntryAccess(matrix_, row, col, ConvertToRef<value_type>());
  }
 
  template <class RowLocalView, class ColLocalView>
  void scatter (const RowLocalView& rowLocalView, const ColLocalView& colLocalView,
                const LocalMatrix& localMatrix, const LocalPattern& localPattern)
  {
    auto rowIndices = Dune::Assembler::Impl::localIndices(rowLocalView.tree());
    auto colIndices = Dune::Assembler::Impl::localIndices(colLocalView.tree());
    localPattern.visitEntries(rowIndices, colIndices, [&](const auto& i, const auto&j) {
      Impl::hybridRecursiveMatrixEntryAccess(matrix_, rowLocalView.index(i), colLocalView.index(j), AccumulateValue(localMatrix[i][j]));
    });
  }
 
  void setZero ()
  {
    matrix_ = 0;
  }
 
private:
  Matrix& matrix_;
};
 
 
 
// Specialization for flat BCRSMatrix
template <class T, class A>
class ISTLMatrixPatternBackend<Dune::BCRSMatrix<T, A>>
  : public Impl::MatrixPatternBackendMixIn<ISTLMatrixPatternBackend<Dune::BCRSMatrix<T, A>>>
{
  using Matrix = Dune::BCRSMatrix<T, A>;
  using Base = Impl::MatrixPatternBackendMixIn<ISTLMatrixPatternBackend<Matrix>>;
 
public:
 
  using LocalPattern = typename Base::LocalPattern;
 
  ISTLMatrixPatternBackend (Matrix& matrix)
    : matrix_(matrix)
  {}
 
  template <class RowBasis, class ColBasis>
  void resize (const RowBasis& rowBasis, const ColBasis& colBasis)
  {
    indexSet_.resize(rowBasis.size(), colBasis.size());
  }
 
  template <class RowIndex, class ColIndex>
  void add (const RowIndex& rowIndex, const ColIndex& colIndex)
  {
    indexSet_.add(rowIndex[0], colIndex[0]);
  }
 
  void setupMatrix ()
  {
    indexSet_.exportIdx(matrix_);
    matrix_ = 0;
  }
 
  template <class F>
  void forEachEntry (F&& f) const
  {
    for (auto row : Dune::range(indexSet_.rows()))
    {
      std::visit([&](const auto& columnIndices) {
        for (auto col : columnIndices)
          f(std::size_t(row), std::size_t(col));
      }, indexSet_.columnIndices(row));
    }
  }
 
private:
  Matrix& matrix_;
  Dune::MatrixIndexSet indexSet_;
};
 
 
 
// Specialization for flat MultiTypeBlockMatrix of flat matrix
template <class Row0, class ... Rows>
class ISTLMatrixPatternBackend<MultiTypeBlockMatrix<Row0,Rows...>>
  : public Impl::MatrixPatternBackendMixIn<ISTLMatrixPatternBackend<MultiTypeBlockMatrix<Row0,Rows...>>>
{
  using Matrix = MultiTypeBlockMatrix<Row0,Rows...>;
  using Base = Impl::MatrixPatternBackendMixIn<ISTLMatrixPatternBackend<Matrix>>;
 
  static constexpr auto outerRows_ = Dune::index_constant<Matrix::N()>();
  static constexpr auto outerCols_ = Dune::index_constant<Matrix::M()>();
 
public:
 
  using LocalPattern = typename Base::LocalPattern;
 
  ISTLMatrixPatternBackend (Matrix& matrix)
    : matrix_(matrix)
  {}
 
  template <class RowBasis, class ColBasis>
  void resize (const RowBasis& rowBasis, const ColBasis& colBasis)
  {
    for (auto i : Dune::range(outerRows_))
      for (auto j : Dune::range(outerCols_))
        indices_[i][j].resize(rowBasis.size({i}), colBasis.size({j}));
  }
 
  template <class RowIndex, class ColIndex>
  void add (const RowIndex& rowIndex, const ColIndex& colIndex)
  {
    indices_[rowIndex[0]][colIndex[0]].add(rowIndex[1], colIndex[1]);
  }
 
  void setupMatrix ()
  {
    Hybrid::forEach(Dune::range(outerRows_), [&](auto&& i) {
      Hybrid::forEach(Dune::range(outerCols_), [&](auto&& j) {
        indices_[i][j].exportIdx(matrix_[i][j]);
        matrix_[i][j] = 0;
      });
    });
  }
 
  template <class F>
  void forEachEntry (F&& f) const
  {
    using MultiIndex = typename Dune::Functions::StaticMultiIndex<std::size_t, 2>;
    Hybrid::forEach(Dune::range(std::size_t(outerRows_)), [&](auto&& i) {
      Hybrid::forEach(Dune::range(std::size_t(outerCols_)), [&](auto&& j) {
        for (auto ii : Dune::range(std::size_t(indices_[i][j].rows())))
        {
          std::visit([&](const auto& columnIndices) {
            for (auto jj : columnIndices)
              f(MultiIndex{{i, ii}}, MultiIndex{{j, std::size_t(jj)}});
          }, indices_[i][j].columnIndices(ii));
        }
      });
    });
  }
 
private:
  Matrix& matrix_;
  std::array<std::array<Dune::MatrixIndexSet,outerCols_>,outerRows_> indices_;
};
 
 
 
// Specialization for Dune::Matrix of flat BCRSMatrix
template <class T, class OuterA, class InnerA>
class ISTLMatrixPatternBackend<Dune::Matrix<Dune::BCRSMatrix<T, InnerA>, OuterA>>
  : public Impl::MatrixPatternBackendMixIn<ISTLMatrixPatternBackend<Dune::Matrix<Dune::BCRSMatrix<T, InnerA>, OuterA>>>
{
  using Matrix = Dune::Matrix<Dune::BCRSMatrix<T, InnerA>, OuterA>;
  using Base = Impl::MatrixPatternBackendMixIn<ISTLMatrixPatternBackend<Matrix>>;
 
public:
 
  using LocalPattern = typename Base::LocalPattern;
 
  ISTLMatrixPatternBackend (Matrix& matrix)
    : matrix_(matrix)
  {}
 
  template <class RowBasis, class ColBasis>
  void resize (const RowBasis& rowBasis, const ColBasis& colBasis)
  {
    auto outerRows = rowBasis.size();
    auto outerCols = colBasis.size();
    indices_.resize(outerRows);
    for (auto i : Dune::range(outerRows))
      indices_[i].resize(outerCols);
    for (auto i : Dune::range(outerRows))
      for (auto j : Dune::range(outerCols))
        indices_[i][j].resize(rowBasis.size({i}), colBasis.size({j}));
  }
 
  template <class RowIndex, class ColIndex>
  void add (const RowIndex& rowIndex, const ColIndex& colIndex)
  {
    indices_[rowIndex[0]][colIndex[0]].add(rowIndex[1], colIndex[1]);
  }
 
  void setupMatrix ()
  {
    auto outerRows = indices_.size();
    auto outerCols = indices_[0].size();
    matrix_.setSize(outerRows, outerCols);
    for (auto i : Dune::range(outerRows))
      for (auto j : Dune::range(outerCols))
        indices_[i][j].exportIdx(matrix_[i][j]);
    for (auto i : Dune::range(outerRows))
      for (auto j : Dune::range(outerCols))
        matrix_[i][j] = 0;
  }
 
  template <class F>
  void forEachEntry (F&& f) const
  {
    using MultiIndex = typename Dune::Functions::StaticMultiIndex<std::size_t, 2>;
    Hybrid::forEach(Dune::range(std::size_t(indices_.size())), [&](auto&& i) {
      Hybrid::forEach(Dune::range(std::size_t(indices_[i].size())), [&](auto&& j) {
        for (auto ii : Dune::range(std::size_t(indices_[i][j].rows())))
        {
          std::visit([&](const auto& columnIndices) {
            for (auto jj : columnIndices)
              f(MultiIndex{{i, ii}}, MultiIndex{{j, std::size_t(jj)}});
          }, indices_[i][j].columnIndices(ii));
        }
      });
    });
  }
 
private:
  Matrix& matrix_;
  std::vector<std::vector<Dune::MatrixIndexSet>> indices_;
};
 
 
 
} // end namespace Dune::Assembler
 
#endif // DUNE_ASSEMBLER_BACKENDS_ISTLMATRIXBACKEND_HH