dense_host_vector.hpp

Go to the documentation of this file.

#pragma once
 
#include <uni20/krylov/matrix_free.hpp>
#include <uni20/krylov/vector_space.hpp>
#include <uni20/tensor/reductions.hpp>
 
#include <cmath>
#include <complex>
#include <cstddef>
#include <span>
#include <stdexcept>
#include <utility>
#include <vector>
 
namespace uni20::krylov
{
 
template <typename Scalar> struct DenseHostVector
{
    std::vector<Scalar> values;
};
 
struct DenseHostVectorTransfer
{
    template <RealOrComplex To, RealOrComplex From>
      requires(Complex<To> == Complex<From>)
    void operator()(DenseHostVector<To>& output, DenseHostVector<From> const& input, Precision p) const
    {
      require_precision<To>(p);
      (void)p.bit_count();
      if (output.values.size() != input.values.size())
        throw std::invalid_argument("Krylov transfer requires matching vector dimensions");
      for (std::size_t i = 0; i < input.values.size(); ++i)
        output.values[i] = gmres_detail::convert_scalar<To>(input.values[i], p);
    }
};
 
template <typename Scalar> class DenseHostVectorAlgebra {
  public:
    explicit DenseHostVectorAlgebra(std::size_t dimension) : dimension_(dimension) {}
 
    [[nodiscard]] std::size_t dimension() const noexcept { return dimension_; }
    [[nodiscard]] std::size_t problem_dimension() const noexcept { return dimension_; }
    [[nodiscard]] std::size_t vector_dimension(DenseHostVector<Scalar> const& x) const noexcept
    {
      return x.values.size();
    }
 
    [[nodiscard]] auto prototype_view(DenseHostVector<Scalar> const& x) const -> std::span<Scalar const>
    {
      return x.values;
    }
    auto prototype_view(DenseHostVector<Scalar> const&&) const = delete;
 
    template <class Element, std::size_t Extent>
      requires std::same_as<std::remove_const_t<Element>, Scalar>
    [[nodiscard]] auto prototype_view(std::span<Element, Extent> x) const -> std::span<Scalar const>
    {
      return x;
    }
 
    [[nodiscard]] std::size_t vector_dimension(std::span<Scalar const> x) const noexcept { return x.size(); }
 
    [[nodiscard]] int allocation_count() const noexcept { return allocation_count_; }
    [[nodiscard]] int copy_count() const noexcept { return copy_count_; }
    [[nodiscard]] int axpy_count() const noexcept { return axpy_count_; }
    [[nodiscard]] int scal_count() const noexcept { return scal_count_; }
    [[nodiscard]] int set_zero_count() const noexcept { return set_zero_count_; }
    [[nodiscard]] int inner_product_count() const noexcept { return inner_product_count_; }
    [[nodiscard]] int norm_count() const noexcept { return norm_count_; }
 
    [[nodiscard]] DenseHostVector<Scalar> allocate_like(DenseHostVector<Scalar> const& x)
    {
      return this->allocate_like(this->prototype_view(x));
    }
 
    [[nodiscard]] DenseHostVector<Scalar> allocate_like(std::span<Scalar const> x)
    {
      ++allocation_count_;
      return DenseHostVector<Scalar>{std::vector<Scalar>(x.size())};
    }
 
    void copy(DenseHostVector<Scalar>& dst, DenseHostVector<Scalar> const& src)
    {
      require_same_size(dst, src);
      ++copy_count_;
      dst.values = src.values;
    }
 
    void axpy(DenseHostVector<Scalar>& y, Scalar alpha, DenseHostVector<Scalar> const& x)
    {
      require_same_size(y, x);
      ++axpy_count_;
      for (std::size_t i = 0; i < y.values.size(); ++i)
      {
        y.values[i] += alpha * x.values[i];
      }
    }
 
    void scal(DenseHostVector<Scalar>& x, Scalar alpha)
    {
      ++scal_count_;
      for (auto& value : x.values)
      {
        value *= alpha;
      }
    }
 
    void set_zero(DenseHostVector<Scalar>& x)
    {
      ++set_zero_count_;
      for (auto& value : x.values)
      {
        value = Scalar{};
      }
    }
 
    [[nodiscard]] Scalar inner_product(DenseHostVector<Scalar> const& x, DenseHostVector<Scalar> const& y)
    {
      require_same_size(x, y);
      ++inner_product_count_;
 
      auto x_span = make_span(x);
      auto y_span = make_span(y);
      return uni20::inner_product_host(uni20::linalg::CpuReferenceBackend{}, x_span, y_span);
    }
 
    [[nodiscard]] uni20::make_real_t<Scalar> norm(DenseHostVector<Scalar> const& x)
    {
      ++norm_count_;
 
      auto x_span = make_span(x);
      return uni20::norm_host(uni20::linalg::CpuReferenceBackend{}, x_span);
    }
 
  private:
    static void require_same_size(DenseHostVector<Scalar> const& lhs, DenseHostVector<Scalar> const& rhs)
    {
      if (lhs.values.size() != rhs.values.size())
      {
        throw std::invalid_argument("dense Krylov vectors have different sizes");
      }
    }
 
    [[nodiscard]] static auto make_span(DenseHostVector<Scalar> const& vector)
    {
      using extents_type = stdex::dextents<uni20::index_type, 1>;
      return stdex::mdspan<Scalar const, extents_type>(
          vector.values.data(), extents_type{static_cast<uni20::index_type>(vector.values.size())});
    }
 
    std::size_t dimension_ = 0;
    int allocation_count_ = 0;
    int copy_count_ = 0;
    int axpy_count_ = 0;
    int scal_count_ = 0;
    int set_zero_count_ = 0;
    int inner_product_count_ = 0;
    int norm_count_ = 0;
};
 
template <typename Scalar> class DenseHostVectorOps : public DenseHostVectorAlgebra<Scalar> {
  public:
    DenseHostVectorOps(std::size_t dimension, std::vector<Scalar> matrix)
        : DenseHostVectorAlgebra<Scalar>(dimension), matrix_(std::move(matrix))
    {
      if (matrix_.size() != dimension * dimension)
        throw std::invalid_argument("dense Krylov operator matrix has the wrong size");
    }
 
    [[nodiscard]] int matvec_count() const noexcept { return matvec_count_; }
 
    void matvec(DenseHostVector<Scalar>& y, DenseHostVector<Scalar> const& x)
    {
      if (x.values.size() != this->problem_dimension() || y.values.size() != this->problem_dimension())
      {
        throw std::invalid_argument("dense Krylov matvec vector has the wrong size");
      }
 
      ++matvec_count_;
      for (std::size_t row = 0; row < this->problem_dimension(); ++row)
      {
        Scalar value{};
        for (std::size_t col = 0; col < this->problem_dimension(); ++col)
        {
          value += matrix_[row * this->problem_dimension() + col] * x.values[col];
        }
        y.values[row] = value;
      }
    }
 
  private:
    std::vector<Scalar> matrix_;
    int matvec_count_ = 0;
};
 
static_assert(KrylovMatrixFreeOperator<DenseHostVectorOps<float>, DenseHostVector<float>, float>);
static_assert(KrylovMatrixFreeOperator<DenseHostVectorOps<double>, DenseHostVector<double>, double>);
static_assert(KrylovMatrixFreeOperator<DenseHostVectorOps<uni20::complex<float>>,
                                       DenseHostVector<uni20::complex<float>>, uni20::complex<float>>);
static_assert(KrylovMatrixFreeOperator<DenseHostVectorOps<uni20::complex<double>>,
                                       DenseHostVector<uni20::complex<double>>, uni20::complex<double>>);
 
} // namespace uni20::krylov