InnerProduct.hpp

#pragma once

#include <functional>
#include <type_traits>
#include <utility>
#include <vector>

#include <dune/common/typeutilities.hh>
#include <amdis/common/ForEach.hpp>
#include <amdis/common/StaticSize.hpp>

namespace AMDiS {
namespace Recursive {

template <class>
struct InnerProduct;

template <class In1, class In2, class T, class BinaryOp1, class BinaryOp2>
T innerProduct (In1 const& in1, In2 const& in2, T init, BinaryOp1 op1, BinaryOp2 op2)
{
  return InnerProduct<In1>::impl(in1, in2, std::move(init), op1, op2);
}


template <class In1, class In2, class T>
T innerProduct (In1 const& in1, In2 const& in2, T init)
{
  return InnerProduct<In1>::impl(in1, in2, std::move(init), std::plus<>{}, std::multiplies<>{});
}


template <class>
struct InnerProduct
{
private:
  // dynamic ranges
  template <class In1, class In2, class T, class BinOp1, class BinOp2,
    class = decltype(std::begin(std::declval<In1>())),
    class = decltype(std::end(std::declval<In1>())),
    class = decltype(std::begin(std::declval<In2>()))>
  static T impl2 (Dune::PriorityTag<2>, In1 const& in1, In2 const& in2, T init,
                  BinOp1 op1, BinOp2 op2)
  {
    auto first1 = std::begin(in1);
    auto first2 = std::begin(in2);
    for (; first1 != std::end(in1); ++first1, ++first2)
      init = Recursive::innerProduct(*first1, *first2, std::move(init), op1, op2);
    return init;
  }

  // ranges with static index access
  template <class In1, class In2, class T, class BinOp1, class BinOp2,
    class = decltype(std::declval<In1>()[std::integral_constant<std::size_t,0>{}])>
  static T impl2 (Dune::PriorityTag<1>, In1 const& in1, In2 const& in2, T init,
                  BinOp1 op1, BinOp2 op2)
  {
    static_assert(static_size_v<In1> == static_size_v<In2>);
    Ranges::forIndices<static_size_v<In1>>([&](auto ii) {
      init = Recursive::innerProduct(in1[ii], in2[ii], std::move(init), op1, op2);
    });
    return init;
  }

  // no range
  template <class In1, class In2, class T, class BinOp1, class BinOp2>
  static T impl2 (Dune::PriorityTag<0>, In1 const& in1, In2 const& in2, T init,
                  BinOp1 op1, BinOp2 op2)
  {
    return op1(std::move(init), op2(in1, in2));
  }

public:
  template <class In1, class In2, class T, class BinOp1, class BinOp2>
  static T impl (In1 const& in1, In2 const& in2, T init, BinOp1 op1, BinOp2 op2)
  {
    return impl2(Dune::PriorityTag<5>{}, in1, in2, init, op1, op2);
  }
};

}} // end namespace AMDiS::Recursive