eigenbackend.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// 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_EIGENBACKEND_HH
#define DUNE_ASSEMBLER_BACKENDS_EIGENBACKEND_HH
 
#include <type_traits>
#include <utility>
#include <vector>
#include <variant>
 
#include <dune/common/indices.hh>
#include <dune/common/rangeutilities.hh>
 
#include <dune/assembler/backends/istlmatrixbackend.hh>
#include <dune/assembler/backends/localpattern.hh>
#include <dune/assembler/backends/flatmatrix.hh>
#include <dune/assembler/backends/flatvector.hh>
#include <dune/assembler/utility.hh>
 
#include <Eigen/Dense>
#include <Eigen/Sparse>
 
 
 
namespace Dune::Assembler {
 
 
 
template <class V>
class EigenVectorBackend;
 
 
 
template <class T>
class EigenVectorBackend<typename Eigen::Matrix<T, Eigen::Dynamic, 1>>
{
public:
  using value_type = T;
 
public:
 
  using Vector = typename Eigen::Matrix<T, Eigen::Dynamic, 1>;
  using LocalVector = Dune::Assembler::FlatVector<value_type>;
 
  EigenVectorBackend (Vector& vector)
    : vector_(vector)
  {}
 
  template <class SizeInfo>
  void resize (const SizeInfo& sizeInfo)
  {
    vector_.resize(sizeInfo.size());
  }
 
  template <class MultiIndex>
  const T& operator[] (const MultiIndex& index) const
  {
    return vector_[index];
  }
 
  template <class MultiIndex>
  T& operator[] (const MultiIndex& index)
  {
    return vector_[index];
  }
 
  template <class LocalView>
  void scatter (const LocalView& localView, const LocalVector& localVector)
  {
    for (std::size_t i = 0; i < localView.size(); ++i)
      (*this)[localView.index(i)] += localVector[i];
  }
 
  template <class TT>
  requires Dune::IsNumber<TT>::value
  EigenVectorBackend& operator= (const TT& value)
  {
    vector_.array() = value;
    return *this;
  }
 
  void setZero ()
  {
    vector_.setZero();
  }
 
  const Vector& vector () const
  {
    return *vector_;
  }
 
  Vector& vector ()
  {
    return *vector_;
  }
 
private:
  Vector& vector_;
};
 
 
 
template <class T>
EigenVectorBackend(Eigen::Matrix<T, Eigen::Dynamic, 1>&)
-> EigenVectorBackend<Eigen::Matrix<T, Eigen::Dynamic, 1>>;
 
 
 
template <class M>
class EigenTripletMatrixPatternBackend;
 
 
 
template <class T>
class EigenTripletMatrixPatternBackend<typename Eigen::SparseMatrix<T, Eigen::RowMajor>>
  : public Impl::MatrixPatternBackendMixIn<EigenTripletMatrixPatternBackend<typename Eigen::SparseMatrix<T, Eigen::RowMajor>>>
{
  using M = Eigen::SparseMatrix<T, Eigen::RowMajor>;
  using Base = Impl::MatrixPatternBackendMixIn<EigenTripletMatrixPatternBackend<M>>;
  using Scalar = typename M::Scalar;
 
public:
 
  using Matrix = Eigen::SparseMatrix<T, Eigen::RowMajor>;
  using LocalPattern = Dune::Assembler::LocalPatternSet;
 
  EigenTripletMatrixPatternBackend (Matrix& matrix)
    : matrix_(matrix)
  {}
 
  template <class RowSizeInfo, class ColSizeInfo>
  void resize (const RowSizeInfo& rowSizeInfo, const ColSizeInfo& colSizeInfo)
  {
    rows_ = rowSizeInfo.size();
    cols_ = colSizeInfo.size();
  }
 
  template <class RowIndex, class ColIndex>
  void add (const RowIndex& rowIndex, const ColIndex& colIndex)
  {
    indexTriplets_.emplace_back(rowIndex[0], colIndex[0], 0);
  }
 
  void setupMatrix ()
  {
    matrix_.resize(rows_, cols_);
    matrix_.setFromTriplets(indexTriplets_.begin(), indexTriplets_.end());
  }
 
private:
  std::vector<Eigen::Triplet<Scalar>> indexTriplets_;
  std::size_t rows_;
  std::size_t cols_;
  Matrix& matrix_;
};
 
 
 
template <class M>
class EigenMatrixPatternBackend;
 
 
 
template <class T>
class EigenMatrixPatternBackend<typename Eigen::SparseMatrix<T, Eigen::RowMajor>>
  : public Impl::MatrixPatternBackendMixIn<EigenMatrixPatternBackend<typename Eigen::SparseMatrix<T, Eigen::RowMajor>>>
{
  using M = Eigen::SparseMatrix<T, Eigen::RowMajor>;
  using Base = Impl::MatrixPatternBackendMixIn<EigenMatrixPatternBackend<M>>;
 
  // Hack: Eigen::SparseMatrix::reserve(rowSizes) tries to detect
  // a range by the presence of decltype(rowSizes)::value_type.
  // Thus we have to augment Dune::TransformedRangeView by a suitable
  // typedef.
  template <class BaseRange, class V>
  class WithValueType
    : public BaseRange
  {
  public:
    using value_type = V;
 
    WithValueType (BaseRange&& baseRange)
      : BaseRange(std::move(baseRange))
    {}
  };
 
  template <class V, class BaseRange>
  static auto withValueType (BaseRange&& baseRange)
  {
    return WithValueType<std::decay_t<BaseRange>, V> (std::move(baseRange));
  }
 
public:
 
  using Matrix = Eigen::SparseMatrix<T, Eigen::RowMajor>;
  using LocalPattern = typename Base::LocalPattern;
 
  EigenMatrixPatternBackend (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 ()
  {
    std::size_t rows = indexSet_.rows();
    std::size_t cols = indexSet_.cols();
 
    matrix_.resize(rows, cols);
 
    auto rowSizes = withValueType<std::size_t>(Dune::transformedRangeView(Dune::range(rows), [&](auto row) { return indexSet_.rowsize(row); }));
    matrix_.reserve(rowSizes);
 
    for (auto row : Dune::range(rows))
      std::visit([&](const auto& colIndices) {
        for (auto col : colIndices)
          matrix_.insert(row, col) = 0;
      }, indexSet_.columnIndices(row));
 
    matrix_.finalize();
  }
 
private:
  Matrix& matrix_;
  Dune::MatrixIndexSet indexSet_;
};
 
 
 
 
template <class M, bool useTriplets=false>
class EigenMatrixBackend;
 
 
 
template <class T, bool useTriplets>
class EigenMatrixBackend<typename Eigen::SparseMatrix<T, Eigen::RowMajor>, useTriplets>
{
 
public:
 
  using Matrix = Eigen::SparseMatrix<T, Eigen::RowMajor>;
  using Entry = typename Matrix::Scalar;
  using value_type = Entry;
  using LocalMatrix = Dune::Assembler::FlatMatrix<value_type>;
  using LocalPattern = Dune::Assembler::LocalPatternSet;
 
  EigenMatrixBackend (Matrix& matrix)
    : matrix_(&matrix)
  {}
 
  EigenMatrixBackend (Matrix& matrix, std::bool_constant<useTriplets>)
    : matrix_(&matrix)
  {}
 
  auto patternBuilder ()
  {
    if constexpr (useTriplets)
      return EigenTripletMatrixPatternBackend<Matrix>(*matrix_);
    else
      return EigenMatrixPatternBackend<Matrix>(*matrix_);
  }
 
  template <class RowMultiIndex, class ColMultiIndex>
  const Entry& operator() (const RowMultiIndex& row, const ColMultiIndex& col) const
  {
    return matrix_->coeffRef(row, col);
  }
 
  template <class RowMultiIndex, class ColMultiIndex>
  Entry& operator() (const RowMultiIndex& row, const ColMultiIndex& col)
  {
    return matrix_->coeffRef(row, col);
  }
 
  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) {
      (*this)(rowLocalView.index(i), colLocalView.index(j)) += localMatrix[i][j];
    });
  }
 
  void setZero ()
  {
    for (auto i: Dune::range(matrix_->rows()))
      for (auto it = Matrix::InnerIterator(*matrix_, i); it; ++it)
        it.valueRef() = 0.0;
  }
 
  const Matrix& matrix () const
  {
    return *matrix_;
  }
 
  Matrix& matrix ()
  {
    return *matrix_;
  }
 
  template <class Value>
  void assign (const Value& value)
  {
    auto& data = matrix_->data();
    for (auto k : Dune::range(data.size()))
      data.value(k) = value;
  }
 
protected:
 
  Matrix* matrix_;
};
 
 
 
template <class T>
EigenMatrixBackend(Eigen::SparseMatrix<T, Eigen::RowMajor>&)
-> EigenMatrixBackend<typename Eigen::SparseMatrix<T, Eigen::RowMajor>>;
 
 
 
} // namespace Dune::Assembler
 
 
 
#endif // DUNE_ASSEMBLER_BACKENDS_EIGENBACKEND_HH