experimental/bcrsmatrix.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_ISTL_EXPERIMENTAL_BCRSMATRIX_HH
#define DUNE_ISTL_EXPERIMENTAL_BCRSMATRIX_HH
 
#include "matrix.hh"
#include "sparseindexranges.hh"
 
#include <dune/common/scalarmatrixview.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/overloadset.hh>
#include <dune/common/std/no_unique_address.hh>
#include <dune/common/std/type_traits.hh>
 
#include <memory>
#include <ranges>
 
namespace Dune::Experimental {
template<class B, class I = uint_least32_t, class A = std::allocator<B>>
class BCRSMatrix : public Matrix<B, SparseIndexRanges<I>, A>
{
  using block_alloc_type = typename std::allocator_traits<A>::template rebind_alloc<B>;
  using block_alloc_traits = std::allocator_traits<block_alloc_type>;
  using void_alloc_type = typename std::allocator_traits<A>::template rebind_alloc<void>;
 
  using SequencedPatternBuilder = SequencedSparseIndexRangeBuilder<void_alloc_type>;
  using UnsequencedPatternBuilder = UnsequencedSparseIndexRangeBuilder<void_alloc_type>;
 
  struct AllocDeleter
  {
    template<class T>
    void operator()(T* builder)
    {
      using rebound_alloc_type = typename std::allocator_traits<A>::template rebind_alloc<T>;
      rebound_alloc_type rebound_alloc(alloc_);
      using rebound_alloc_traits = typename std::allocator_traits<rebound_alloc_type>;
      rebound_alloc_traits::destroy(rebound_alloc, builder);
      rebound_alloc_traits::deallocate(rebound_alloc, builder, 1);
    }
    DUNE_NO_UNIQUE_ADDRESS A alloc_;
  };
 
  using BuilderVariant = std::variant<
    std::monostate,
    std::unique_ptr<SequencedPatternBuilder, AllocDeleter>,
    std::unique_ptr<UnsequencedPatternBuilder, AllocDeleter>>;
 
  class DefaultBlockAssgnment
  {
    DefaultBlockAssgnment() = default;
    friend class BCRSMatrix;
  public:
    DefaultBlockAssgnment& operator=(const B& e);
    DefaultBlockAssgnment& operator=(const typename FieldTraits<B>::field_type& f)
    requires (!std::same_as<B, typename FieldTraits<B>::field_type>);
  };
 
public:
  using typename Matrix<B, SparseIndexRanges<I>, A>::block_type;
  using typename Matrix<B, SparseIndexRanges<I>, A>::size_type;
  using typename Matrix<B, SparseIndexRanges<I>, A>::pattern_type;
  using allocator_type = A;
 
  enum BuildMode : uint_least8_t
  {
    row_wise = 0,
    random = 1,
    implicit = 1,
    unknown
  };
 
  struct CompressionStatistics
  {
    double avg;
    size_type maximum;
    size_type overflow_total;
    double mem_ratio;
  };
 
  // inherit constructors TODO check that this works
  using Matrix<B, SparseIndexRanges<I>, A>::Matrix;
 
  using Matrix<B, SparseIndexRanges<I>, A>::operator=;
 
  using Matrix<B, SparseIndexRanges<I>, A>::get_allocator;
 
  BCRSMatrix(const Matrix<B, SparseIndexRanges<I>, A>& other,
              const allocator_type& allocator = allocator_type())
    : Matrix<B, SparseIndexRanges<I>, A>(other, allocator)
  {
  }
 
  BCRSMatrix(Matrix<B, SparseIndexRanges<I>, A>&& other,
              const allocator_type& allocator = allocator_type())
    : Matrix<B, SparseIndexRanges<I>, A>(std::move(other), allocator)
  {
  }
 
  BCRSMatrix(const BCRSMatrix& other,
              const allocator_type& allocator = allocator_type());
  BCRSMatrix(BCRSMatrix&& other,
              const allocator_type& allocator = allocator_type());
 
  BCRSMatrix<B, I, A>& operator=(const BCRSMatrix& other);
 
  BCRSMatrix<B, I, A>& operator=(BCRSMatrix&& other);
 
  BCRSMatrix(size_type rows,
             size_type cols,
             BuildMode bm,
             const allocator_type& allocator = allocator_type());
 
  BCRSMatrix(size_type rows,
             size_type cols,
             size_type nnz,
             BuildMode bm,
             const allocator_type& allocator = allocator_type());
 
  BCRSMatrix(size_type rows,
             size_type cols,
             size_type avg,
             double overflowRatio,
             BuildMode bm,
             const allocator_type& allocator = allocator_type());
 
  void swap(BCRSMatrix& other)
  {
    using std::swap;
    swap(static_cast<Matrix<B, SparseIndexRanges<I>, A>&>(*this),
         static_cast<Matrix<B, SparseIndexRanges<I>, A>&>(other));
    swap(builder_, other.builder_);
  }
 
  friend void swap(BCRSMatrix& lhs, BCRSMatrix& rhs)
  {
    lhs.swap(rhs);
  }
 
  void setBuildMode(BuildMode bm);
 
  BuildMode buildMode() const;
 
  void setSize(size_type rows, size_type cols, size_type nnz = 0);
 
 
  void setrowsize(size_type row, size_type size);
 
  void incrementrowsize(size_type row, size_type increment = 1);
 
  void endrowsizes() {}
 
  void addindex(size_type row, size_type col);
  void endindices();
 
  template<typename It>
  void setIndicesNoSort(size_type row, It begin, It end) {
    DUNE_THROW(NotImplemented, "setIndicesNoSort is not implemented yet for BCRSMatrix row-wise build mode");
  }
 
  template<typename It>
  void setIndices(size_type row, It begin, It end) {
    DUNE_THROW(NotImplemented, "setIndices is not implemented yet for BCRSMatrix row-wise build mode");
  }
 
 
  class CreateIterator
    : public SequencedPatternBuilder::IndexRangeBuilder
  {
    using Base = typename SequencedPatternBuilder::IndexRangeBuilder;
    using EndSentinel = typename SequencedPatternBuilder::EndSentinel;
 
    CreateIterator(Base&& pattern_builder, BCRSMatrix* sm);
 
    CreateIterator(const CreateIterator&) = delete;
    CreateIterator& operator=(const CreateIterator&) = delete;
  public:
    CreateIterator(CreateIterator&& other) noexcept;
    CreateIterator& operator=(CreateIterator&& other) noexcept;
    CreateIterator& operator++();
 
  private:
    friend class BCRSMatrix;
    BCRSMatrix* matrix_;
  };
 
  CreateIterator createbegin();
 
  typename SequencedPatternBuilder::EndSentinel createend();
 
  DefaultBlockAssgnment entry(size_type row, size_type col);
 
  CompressionStatistics compress();
 
private:
  template<class Builder>
  static Builder& getBuilder(BuilderVariant& builder,
                             std::string_view build_mode_name)
  {
    using BuilderPtrType = std::unique_ptr<Builder, AllocDeleter>;
    if (!std::holds_alternative<BuilderPtrType>(builder))
      DUNE_THROW(
        InvalidStateException,
        "Builder "
          << className<Builder>()
          << " can only be used when building a BCRSMatrix with BuildMode::"
          << build_mode_name);
 
    const auto& ptr = std::get<BuilderPtrType>(builder);
    if (!ptr)
      DUNE_THROW(
        InvalidStateException,
        "Builder "
          << className<Builder>()
          << " can only be called once the builder has been initialized");
    return *ptr;
  }
 
  SequencedPatternBuilder& getSequencedPatternBuilder()
  {
    return getBuilder<SequencedPatternBuilder>(builder_, "row_wise");
  }
 
  UnsequencedPatternBuilder& getUnsequencedPatternBuilder()
  {
    return getBuilder<UnsequencedPatternBuilder>(builder_, "implicit");
  }
 
  BuilderVariant builder_;
};
 
// Matrix implementation
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>::BCRSMatrix(size_type n,
                                  size_type m,
                                  size_type nnz,
                                  BuildMode bm,
                                  const allocator_type& allocator)
  : Matrix<B, SparseIndexRanges<I>, A>(nullptr, allocator)
{
  setBuildMode(bm);
  setSize(n, m, nnz);
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>::BCRSMatrix(size_type n,
                                  size_type m,
                                  BuildMode bm,
                                  const allocator_type& allocator)
  : BCRSMatrix(n, m, 0, bm, allocator)
{
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>::BCRSMatrix(size_type n,
                                  size_type m,
                                  size_type avg,
                                  double overflowRatio,
                                  BuildMode bm,
                                  const allocator_type& allocator)
  : BCRSMatrix(n,
                m,
                static_cast<size_type>((avg * n) * overflowRatio),
                bm,
                allocator)
{
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>::BCRSMatrix(const BCRSMatrix& other,
                                  const allocator_type& allocator)
  : Matrix<B, SparseIndexRanges<I>, A>(other, allocator)
{
  builder_ = std::monostate{};
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>::BCRSMatrix(BCRSMatrix&& other,
                                  const allocator_type& allocator)
  : Matrix<B, SparseIndexRanges<I>, A>(std::move(other), allocator)
{
  builder_ = std::move(other.builder_);
}
 
template<class B, class I, class A>
void
BCRSMatrix<B, I, A>::setBuildMode(BuildMode bm)
{
  *this = BCRSMatrix(get_allocator());
  switch (bm) {
    case BuildMode::row_wise:
      builder_ =
        std::unique_ptr<SequencedPatternBuilder, AllocDeleter>(
          nullptr, AllocDeleter{ get_allocator() });
      break;
    case BuildMode::implicit: // same as BuildMode::random
      builder_ =
        std::unique_ptr<UnsequencedPatternBuilder, AllocDeleter>(
          nullptr, AllocDeleter{ get_allocator() });
      break;
    case BuildMode::unknown:
      builder_ = std::monostate{};
      break;
  }
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::buildMode() const -> BuildMode
{
  return std::visit(
    overload(
      [](std::monostate) { return BuildMode::unknown; },
      [](const std::unique_ptr<SequencedPatternBuilder,
                               AllocDeleter>&) { return BuildMode::row_wise; },
      [](const std::unique_ptr<UnsequencedPatternBuilder,
                               AllocDeleter>&) { return BuildMode::implicit; }),
    builder_);
}
 
template<class B, class I, class A>
void
BCRSMatrix<B, I, A>::setSize(size_type rows, size_type cols, size_type nnz)
{
  using seq_builder_alloc_type = typename std::allocator_traits<
    allocator_type>::template rebind_alloc<SequencedPatternBuilder>;
  using unseq_builder_alloc_type = typename std::allocator_traits<
    allocator_type>::template rebind_alloc<UnsequencedPatternBuilder>;
  const auto void_alloc = void_alloc_type(get_allocator());
  size_type avg = (rows > 0) ? (nnz / rows) : 0;
  std::visit(
    overload(
      [](std::monostate) {
        DUNE_THROW(InvalidStateException,
                   "Cannot set size of BCRSMatrix with unknown build mode");
      },
      [balloc = seq_builder_alloc_type(get_allocator()),
       void_alloc,
       rows,
       cols,
       avg](std::unique_ptr<SequencedPatternBuilder, AllocDeleter>&
              builder) mutable {
        builder.reset(
          std::allocator_traits<seq_builder_alloc_type>::allocate(balloc, 1));
        std::allocator_traits<seq_builder_alloc_type>::construct(
          balloc, builder.get(), rows, cols, avg, void_alloc);
      },
      [balloc = unseq_builder_alloc_type(get_allocator()),
       void_alloc,
       rows,
       cols,
       avg](std::unique_ptr<UnsequencedPatternBuilder, AllocDeleter>&
              builder) mutable {
        builder.reset(
          std::allocator_traits<unseq_builder_alloc_type>::allocate(balloc, 1));
        std::allocator_traits<unseq_builder_alloc_type>::construct(
          balloc, builder.get(), rows, cols, avg, void_alloc);
      }),
    builder_);
}
 
template<class B, class I, class A>
void
BCRSMatrix<B, I, A>::setrowsize(size_type row, size_type size)
{
  getUnsequencedPatternBuilder().setSize(row, size);
}
 
template<class B, class I, class A>
void
BCRSMatrix<B, I, A>::incrementrowsize(size_type row, size_type increment)
{
  getUnsequencedPatternBuilder().incrementSize(row, increment);
}
 
template<class B, class I, class A>
void
BCRSMatrix<B, I, A>::addindex(size_type row, size_type col)
{
  getUnsequencedPatternBuilder().addIndex(row, col);
}
 
template<class B, class I, class A>
void
BCRSMatrix<B, I, A>::endindices()
{
  using pattern_alloc_type = typename std::allocator_traits<
    allocator_type>::template rebind_alloc<SparseIndexRanges<I>>;
  pattern_alloc_type pattern_alloc(get_allocator());
  std::shared_ptr<pattern_type> pt = std::allocate_shared<pattern_type>(
    pattern_alloc, std::move(getUnsequencedPatternBuilder()));
  *this = Matrix<B, SparseIndexRanges<I>, A>(std::move(pt), get_allocator());
  builder_ = std::monostate{};
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>::CreateIterator::CreateIterator(Base&& pattern_builder,
                                                     BCRSMatrix* sm)
  : Base(std::move(pattern_builder))
  , matrix_(sm)
{
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>::CreateIterator::CreateIterator(
  CreateIterator&& other) noexcept
  : Base(std::move(other))
  , matrix_(std::exchange(other.matrix_, nullptr))
{
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::CreateIterator::operator=(CreateIterator&& other) noexcept
  -> CreateIterator&
{
  if (this != &other) {
    Base::operator=(std::move(other));
    matrix_ = std::exchange(other.matrix_, nullptr);
  }
  return *this;
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::CreateIterator::operator++() -> CreateIterator&
{
  Base::operator++();
  assert(matrix_);
  if (*this == EndSentinel()) {
    // build pattern and allocate data when we have reached the end of the index
    // sets
    using pattern_alloc_type =
      typename block_alloc_traits::template rebind_alloc<pattern_type>;
    pattern_alloc_type pattern_alloc(matrix_->get_allocator());
    std::shared_ptr<pattern_type> pt = std::allocate_shared<pattern_type>(
      pattern_alloc, std::move(matrix_->getSequencedPatternBuilder()));
    *matrix_ =
      Matrix<B, SparseIndexRanges<I>, A>(std::move(pt), matrix_->get_allocator());
    matrix_->builder_ = std::monostate{};
  }
  return *this;
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::createbegin() -> CreateIterator
{
  return { getSequencedPatternBuilder().begin(), this };
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::createend() ->
  typename SequencedPatternBuilder::EndSentinel
{
  return {};
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::entry(size_type row, size_type col)
  -> DefaultBlockAssgnment
{
  getUnsequencedPatternBuilder().addIndex(row, col);
  return {};
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::compress() -> CompressionStatistics
{
  endindices();
 
  size_type max_pattern_size = 0;
  for (const auto& index_range : this->pattern())
    max_pattern_size =
      std::max<size_type>(max_pattern_size, index_range.size());
  return CompressionStatistics{ /*avg = */ this->pattern().count() /
                                  static_cast<double>(this->pattern().size()),
                                /* maximum = */ max_pattern_size,
                                /* overflow_total = */ 0,
                                /* mem_ratio = */ 1.0 };
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::DefaultBlockAssgnment::operator=(const B& e)
  -> DefaultBlockAssgnment&
{
  if (e != B{})
    DUNE_THROW(RangeError,
               "BCRSMatrix with BuildMode::row_wise only supports default "
               "constructed entries, but got "
                 << e);
  return *this;
}
 
template<class B, class I, class A>
auto
BCRSMatrix<B, I, A>::DefaultBlockAssgnment::operator=(const typename FieldTraits<B>::field_type& e)
-> DefaultBlockAssgnment&
requires (!std::same_as<B, typename FieldTraits<B>::field_type>)
{
  if (e != typename FieldTraits<B>::field_type{})
    DUNE_THROW(RangeError,
               "BCRSMatrix with BuildMode::row_wise only supports default "
               "constructed entries, but got "
                 << e);
  return *this;
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>&
BCRSMatrix<B, I, A>::operator=(const BCRSMatrix& other)
{
  if (this != &other) {
    static_cast<Matrix<B, SparseIndexRanges<I>, A>&>(*this) =
      static_cast<const Matrix<B, SparseIndexRanges<I>, A>&>(other);
    builder_ = std::monostate{};
  }
  return *this;
}
 
template<class B, class I, class A>
BCRSMatrix<B, I, A>&
BCRSMatrix<B, I, A>::operator=(BCRSMatrix&& other)
{
  if (this != &other) {
    static_cast<Matrix<B, SparseIndexRanges<I>, A>&>(*this) = std::move(other);
    builder_ = std::move(other.builder_);
  }
  return *this;
}
 
} //namespace Dune::Experimental
 
namespace Dune {
template<class B, class I, class A>
struct FieldTraits<Experimental::BCRSMatrix<B, I, A>>
  : public FieldTraits<typename Experimental::BCRSMatrix<B, I, A>::field_type>
{};
 
} // namespace Dune
 
#endif // DUNE_ISTL_EXPERIMENTAL_BCRSMATRIX_HH