19#ifndef UTILS_MATRIX_VECTOR_PRODUCT_HPP
20#define UTILS_MATRIX_VECTOR_PRODUCT_HPP
32template <
int c,
typename T>
struct mul {
33 constexpr T operator()(
const T a)
const {
return c * a; }
36template <
typename T>
struct mul<0, T> {
37 constexpr T operator()(
const T)
const {
return T{}; }
40template <
typename T>
struct mul<1, T> {
41 constexpr T operator()(
const T a)
const {
return a; }
44template <
int I,
typename T,
typename Container, std::size_t N,
int c,
46struct inner_product_template_impl {
47 constexpr T operator()(
const Container &vec)
const {
48 return mul<c, T>{}(get<I>(vec)) +
49 inner_product_template_impl<I + 1, T, Container, N, cs...>{}(vec);
53template <
int I,
typename T,
typename Container, std::
size_t N,
int c>
54struct inner_product_template_impl<I, T, Container, N, c> {
55 constexpr T operator()(
const Container &vec)
const {
56 return mul<c, T>{}(get<I>(vec));
60template <
typename T, std::size_t N,
61 const std::array<std::array<int, N>, N> &matrix,
typename Container,
62 std::size_t row_index, std::size_t... column_indices>
63constexpr T inner_product_helper(
const Container &vec,
64 std::index_sequence<column_indices...>) {
65 return inner_product_template_impl<
66 0, T, Container, N, get<column_indices>(get<row_index>(matrix))...>{}(
70template <
typename T, std::size_t N,
71 const std::array<std::array<int, N>, N> &matrix,
typename Container,
72 std::size_t row_index>
73constexpr T inner_product_template(
const Container &vec) {
74 return detail::inner_product_helper<T, N, matrix, Container, row_index>(
75 vec, std::make_index_sequence<N>{});
78template <
typename T, std::size_t N,
79 const std::array<std::array<int, N>, N> &matrix,
typename Container,
80 std::size_t... column_indices>
81constexpr std::array<T, N>
82matrix_vector_product_helper(
const Container &vec,
83 std::index_sequence<column_indices...>) {
84 return std::array<T, N>{
85 {inner_product_template<T, N, matrix, Container, column_indices>(
102template <
typename T, std::size_t N,
103 const std::array<std::array<int, N>, N> &matrix,
typename Container>
105 return detail::matrix_vector_product_helper<T, N, matrix, Container>(
106 vec, std::make_index_sequence<N>{});
constexpr std::array< T, N > matrix_vector_product(const Container &vec)
Calculate the matrix-vector product for a statically given (square) matrix.