backends/cpu/matrix_exponential.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <uni20/config.hpp>
 
#if UNI20_HAS_FLOAT128 && UNI20_FLOAT128_PROVIDER_MPLAPACK
#include <mplapack_config.h>
#endif
 
#include <uni20/common/trace.hpp>
#include <uni20/core/scalar_concepts.hpp>
#include <uni20/linalg/dispatch.hpp>
#include <uni20/linalg/operation_tags.hpp>
#include <uni20/mdspan/concepts.hpp>
#include <uni20/tensor/access.hpp>
#include <uni20/tensor/concepts.hpp>
#include <uni20/tensor/tensor.hpp>
 
#include <complex>
#include <concepts>
#include <type_traits>
 
namespace uni20::linalg::backends::cpu
{
 
template <uni20::RealOrComplex Scalar>
uni20::DenseMatrix<Scalar> matrix_exponential(uni20::DenseMatrix<Scalar> const& matrix, uni20::make_real_t<Scalar> t);
 
template <uni20::RealOrComplex Scalar>
uni20::DenseMatrix<uni20::complex<uni20::make_real_t<Scalar>>>
matrix_exponential(uni20::DenseMatrix<Scalar> const& matrix, uni20::complex<uni20::make_real_t<Scalar>> t);
 
} // namespace uni20::linalg::backends::cpu
 
namespace uni20::linalg
{
namespace detail::cpu_reference
{
 
template <uni20::MutableRankedMdspanLike<2> OutputMdspan, uni20::RankedMdspanLike<2> InputMdspan, class TimeScalar>
consteval auto matrix_exponential_acceptance()
{
  using input_scalar = std::remove_cv_t<typename InputMdspan::element_type>;
  using input_matrix = uni20::DenseMatrix<input_scalar>;
 
  if constexpr (requires(input_matrix const& input, TimeScalar time) {
                  backends::cpu::matrix_exponential(input, time);
                })
  {
    using result_matrix =
        decltype(backends::cpu::matrix_exponential(std::declval<input_matrix const&>(), std::declval<TimeScalar>()));
    if constexpr (requires(result_matrix& result, OutputMdspan& output, InputMdspan& input,
                           typename OutputMdspan::index_type output_index,
                           typename InputMdspan::index_type input_index) {
                    output.operator[](output_index, output_index) = result[std::size_t{}, std::size_t{}];
                    static_cast<input_scalar>(input.operator[](input_index, input_index));
                  })
    {
      return kernel_types_yes;
    }
    else
    {
      return kernel_types_no;
    }
  }
  else
  {
    return kernel_types_no;
  }
}
 
template <class OutputMdspan, class InputMdspan, class TimeScalar>
KernelAttempt matrix_exponential(OutputMdspan& output, InputMdspan& input, TimeScalar time)
{
  using input_type = std::remove_cvref_t<InputMdspan>;
  using input_scalar = std::remove_cv_t<typename input_type::element_type>;
  using input_index_type = typename input_type::index_type;
  using output_index_type = typename std::remove_cvref_t<OutputMdspan>::index_type;
 
  CHECK_EQUAL(input.extent(0), input.extent(1));
  CHECK_EQUAL(output.extent(0), input.extent(0));
  CHECK_EQUAL(output.extent(1), input.extent(1));
 
  std::size_t const rows = static_cast<std::size_t>(input.extent(0));
  std::size_t const cols = static_cast<std::size_t>(input.extent(1));
  uni20::DenseMatrix<input_scalar> materialized(uni20::uninitialized, rows, cols);
  for (input_index_type row = 0; row < input.extent(0); ++row)
  {
    for (input_index_type col = 0; col < input.extent(1); ++col)
    {
      materialized[static_cast<std::size_t>(row), static_cast<std::size_t>(col)] =
          static_cast<input_scalar>(input[row, col]);
    }
  }
 
  auto result = backends::cpu::matrix_exponential(materialized, time);
  for (output_index_type row = 0; row < output.extent(0); ++row)
  {
    for (output_index_type col = 0; col < output.extent(1); ++col)
    {
      output[row, col] = result[static_cast<std::size_t>(row), static_cast<std::size_t>(col)];
    }
  }
  return KernelAttempt::success;
}
} // namespace detail::cpu_reference
 
template <uni20::MutableRankedMdspecLike<2> OutputMdspan, uni20::RankedMdspecLike<2> InputMdspan, class TimeScalar>
  requires uni20::HostWritableMdspec<OutputMdspan> && uni20::HostReadableMdspec<InputMdspan>
consteval auto kernel_accepts_types(CpuReferenceBackend const&, matrix_exponential_op const&, OutputMdspan&,
                                    InputMdspan&, TimeScalar const&)
{
  using output_span = uni20::host_write_mdspan_t<OutputMdspan>;
  using input_span = uni20::host_read_mdspan_t<InputMdspan>;
  constexpr auto acceptance =
      detail::cpu_reference::matrix_exponential_acceptance<output_span, input_span, TimeScalar>();
  if constexpr (acceptance == KernelTypeAcceptance::yes)
    return kernel_types_yes;
  else
    return kernel_types_no;
}
 
template <uni20::MutableRankedMdspecLike<2> OutputMdspan, uni20::RankedMdspecLike<2> InputMdspan, class TimeScalar>
  requires uni20::HostWritableMdspec<OutputMdspan> && uni20::HostReadableMdspec<InputMdspan>
KernelAttempt try_kernel(CpuReferenceBackend, matrix_exponential_op const&, OutputMdspan& output, InputMdspan& input,
                         TimeScalar time)
{
  auto output_access = acquire_host_write_access_sync(output);
  auto input_access = acquire_host_read_access_sync(input);
  auto output_span = output_access.mdspan();
  auto input_span = input_access.mdspan();
  return detail::cpu_reference::matrix_exponential(output_span, input_span, time);
}
 
} // namespace uni20::linalg