InnerProduct.hpp 3.14 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#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;

18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
/**
 * \brief Recursive inner product of two containers [in1] and [in2] following the
 * signature of `std::inner_product` but applied to the leaf elements of the possibly
 * hierarchic containers.
 *
 * This implements the reduction operation of the form
 * ```
 * T result = init;
 * for (i...)
 *   result = op1(result, op2(in1[i], in2[i]);
 * return result;
 * ```
 **/
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)
33
34
35
36
{
  return InnerProduct<In1>::impl(in1, in2, std::move(init), op1, op2);
}

37
38
39
40
41
42
43
44
45
46
47
/// Specialization of \ref innerProduct if no functors are given, use plus and multiplies
/// by default.
/**
 * This implements the standard Euclidean inner product
 * ```
 * result = init + sum_i in1[i] * in2[i]
 * ```
 *
 * Note, no conjugation of the first or second argument is used in this default.
 **/

48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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<>{});
}


/// General implementation of recursive inner-product
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>()))>
65
66
  static T impl2 (Dune::PriorityTag<2>, In1 const& in1, In2 const& in2, T init,
                  BinOp1 op1, BinOp2 op2)
67
68
69
70
71
72
73
74
75
76
77
  {
    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>{}])>
78
79
  static T impl2 (Dune::PriorityTag<1>, In1 const& in1, In2 const& in2, T init,
                  BinOp1 op1, BinOp2 op2)
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
  {
    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