numpybackend.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-FUFEM 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_FUFEM_BACKENDS_NUMPYBACKEND_HH
#define DUNE_FUFEM_BACKENDS_NUMPYBACKEND_HH
 
#include <cstddef>
#include <type_traits>
#include <utility>
#include <algorithm>
#include <array>
#include <tuple>
 
#include <dune/common/rangeutilities.hh>
 
#include <dune/assembler/backends/flatvector.hh>
#include <dune/assembler/backends/flatmatrix.hh>
#include <dune/assembler/backends/localpattern.hh>
#include <dune/assembler/utility.hh>
 
#define DUNE_FUFEM_DISABLE_HEADER_DEPRECATION
#include <dune/fufem/backends/matrixbuilder.hh>
#undef DUNE_FUFEM_DISABLE_HEADER_DEPRECATION
 
#include <dune/python/pybind11/numpy.h>
 
 
 
namespace Dune::Fufem {
 
template<class T>
class NumPyVectorBackend
{
 
public:
 
  using Vector = pybind11::array_t<T>;
 
  NumPyVectorBackend(Vector& vector) :
    vector_(&vector)
  {}
 
  template<class SizeInfo>
  void resize(const SizeInfo& sizeInfo)
  {
    *vector_ = pybind11::array_t<T>(sizeInfo.size());
  }
 
  template<class MultiIndex>
  decltype(auto) operator[](const MultiIndex& index) const
  {
    auto vector_u = vector_->unchecked();
    return vector_u[index];
  }
 
  template<class MultiIndex>
  decltype(auto) operator[](const MultiIndex& index)
  {
    auto vector_mu = vector_->mutable_unchecked();
    return vector_mu[index];
  }
 
  template<typename Value>
  void operator= (const Value& value)
  {
    auto size = vector_->shape(0);
    auto vector_mu = vector_->mutable_unchecked();
    for(auto&& i : Dune::range(size))
      vector_mu(i) = value;
  }
 
  void operator= (const NumPyVectorBackend<T>& other)
  {
    vector() = other.vector();
  }
 
  const Vector& vector() const
  {
    return *vector_;
  }
 
  Vector& vector()
  {
    return *vector_;
  }
 
  // Extended interface of a Dune::Assembler vector-backend
 
  using value_type = T;
  using LocalVector = Dune::Assembler::FlatVector<value_type>;
 
  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];
  }
 
  void setZero ()
  {
    (*this) = 0;
  }
 
private:
  Vector* vector_;
};
 
 
 
template<class T=double>
struct NumPyCSRMatrix
{
  using size_type = std::size_t;
  using field_type = T;
 
  pybind11::array_t<T> entries;
  pybind11::array_t<size_type> colIndices;
  pybind11::array_t<size_type> rowPtrs;
 
  std::array<size_type, 2> shape;
 
  auto asTuple()
  {
    return std::tuple(std::tuple(entries, colIndices, rowPtrs), shape);
  }
 
  auto asTuple() const
  {
    return std::tuple(std::tuple(entries, colIndices, rowPtrs), shape);
  }
 
};
 
 
 
namespace Impl {
 
  // This implements the pattern builder for NumPyCSRMatrix<T>
  template<class T>
  class NumPyCSRMatrixPatternBuilder
  {
  public:
 
    using Matrix = NumPyCSRMatrix<T>;
    using LocalPattern = Dune::Assembler::LocalPatternMatrix;
    using size_type = typename Matrix::size_type;
 
    NumPyCSRMatrixPatternBuilder(Matrix& matrix) :
      matrix_(matrix)
    {}
 
    template<class RowSizeInfo, class ColSizeInfo>
    void resize(const RowSizeInfo& rowSizeInfo, const ColSizeInfo& colSizeInfo)
    {
      indices_.resize(rowSizeInfo.size(), colSizeInfo.size());
    }
 
    void setupMatrix()
    {
      size_type rows = indices_.rows();
      size_type cols = indices_.cols();
      matrix_.shape = {rows, cols};
 
      // Compute number of nonzeros
      std::size_t nnz = 0;
      for(auto row : Dune::range(rows))
        nnz += indices_.rowsize(row);
 
      // Allocate data
      matrix_.entries = pybind11::array_t<T>(nnz);
      matrix_.colIndices = pybind11::array_t<size_type>(nnz);
      matrix_.rowPtrs = pybind11::array_t<size_type>(rows+1);
 
      auto entries_mu = matrix_.entries.mutable_unchecked();
      auto colIndices_mu = matrix_.colIndices.mutable_unchecked();
      auto rowPtrs_mu = matrix_.rowPtrs.mutable_unchecked();
 
      // Write column indices
      size_type next = 0;
      for(auto row : Dune::range(rows))
      {
        rowPtrs_mu(row) = next;
        std::visit([&](const auto& colIndicesOfRow) {
          for(auto&& colIndex : colIndicesOfRow)
          {
            colIndices_mu(next) = colIndex;
            ++next;
          }
        }, indices_.columnIndices(row));
      }
      rowPtrs_mu(rows) = next;
 
      // Zero initialize entries
      for(auto i : Dune::range(nnz))
        entries_mu(i) = 0;
    }
 
    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) {
        add(rowLocalView.index(i), colLocalView.index(j));
      });
    }
 
    template <class F>
    void forEachEntry (F&& f) const
    {
      for (auto row : Dune::range(indices_.rows()))
      {
        std::visit([&](const auto& columnIndices) {
          for (auto col : columnIndices)
            f(std::size_t(row), std::size_t(col));
        }, indices_.columnIndices(row));
      }
    }
 
    // Old dune-fufem interface
 
    template <class RowIndex, class ColIndex>
    void add (const RowIndex& rowIndex, const ColIndex& colIndex)
    {
      indices_.add(rowIndex[0], colIndex[0]);
    }
 
    template<class RowIndex, class ColIndex>
    void insertEntry(const RowIndex& rowIndex, const ColIndex& colIndex)
    {
      this->add(rowIndex, colIndex);
    }
 
  private:
    Dune::MatrixIndexSet indices_;
    Matrix& matrix_;
  };
 
}
 
 
 
 
template<class T=double>
class NumPyCSRMatrixBackend
{
public:
 
  using value_type = T;
  using LocalMatrix = Dune::Assembler::FlatMatrix<value_type>;
  using PatternBuilder = Impl::NumPyCSRMatrixPatternBuilder<T>;
  using LocalPattern = typename PatternBuilder::LocalPattern;
 
  using Matrix = NumPyCSRMatrix<T>;
 
  NumPyCSRMatrixBackend(Matrix& matrix) :
    matrix_(&matrix)
  {}
 
  PatternBuilder patternBuilder()
  {
    return {*matrix_};
  }
 
  template<class RowIndex, class ColIndex>
  auto nnzIndex(const RowIndex& row, const ColIndex& col) const
  {
    auto colIndices_u = matrix_->colIndices.unchecked();
    auto rowPtrs_u = matrix_->rowPtrs.unchecked();
    auto ptrRange = Dune::range(rowPtrs_u(row), rowPtrs_u(row+1));
    auto colIndicesOfRow = Dune::transformedRangeView(ptrRange, [&](auto i) -> decltype(auto) { return colIndices_u(i); });
    auto k = std::lower_bound(colIndicesOfRow.begin(), colIndicesOfRow.end(), col) - colIndicesOfRow.begin();
    return k+rowPtrs_u(row);
  }
 
  template<class RowIndex, class ColIndex>
  const value_type& operator()(const RowIndex& row, const ColIndex& col) const
  {
    auto entries_u = matrix_->entries.unchecked();
    return entries_u(nnzIndex(row, col));
  }
 
  template<class RowIndex, class ColIndex>
  value_type& operator()(const RowIndex& row, const ColIndex& col)
  {
    auto entries_mu = matrix_->entries.mutable_unchecked();
    return entries_mu(nnzIndex(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 ()
  {
    assign(0);
  }
 
  // Old dune-fufem interface
 
  using Entry = T;
 
  const Matrix& matrix() const
  {
    return *matrix_;
  }
 
  Matrix& matrix()
  {
    return *matrix_;
  }
 
  template<class Value>
  void assign(const Value& value)
  {
    auto nnz = matrix_->entries.shape(0);
    auto entries_mu = matrix_->entries.mutable_unchecked();
    for(auto&& i : Dune::range(nnz))
      entries_mu(i) = value;
  }
 
protected:
 
  Matrix* matrix_;
};
 
 
 
template<class T>
class
[[deprecated("This class is deprecated and will be removed after 2.11.")]]
MatrixBuilder<NumPyCSRMatrix<T>>
 : public Impl::NumPyCSRMatrixPatternBuilder<T>
{};
 
} // namespace Dune::Fufem
 
 
 
#endif // DUNE_FUFEM_BACKENDS_NUMPYBACKEND_HH