D3732R0: Parallel numeric ranges algorithms: Performance and parallelism

[mhoemmen]

D3732R0: Parallel numeric ranges algorithms: Performance and parallelism

Executive summary

What we propose

Ranges overloads (both parallel and nonparallel) of the following algorithms:

• reduce, and unary and binary transform_reduce
• inclusive_scan and transform_inclusive_scan
• exclusive_scan and transform_exclusive_scan

Parallel and non-parallel convenience wrappers:

• ranges::sum and ranges::product for special cases of reduce with addition and multiplication, respectively; and
• ranges::dot for the special case of binary transform_reduce with transform multiplies{} and reduction plus{}

With the following features:

• Ranges as output (like P3179)
• Parallel algorithms take sized random access ranges
• Non-parallel algorithms take sized forward ranges

Performance issues affecting design

1. Lack of an identity value complicates parallelization
2. *transform_view not trivially copyable even if function is
3. (Ranges) views generally can hinder optimizations
4. Expose built-in arithmetic on arithmetic types, where possible

Questions for SG1

1. Do you agree that we need transform_* versions of algorithms?
2. Do you agree that we need a way for users to specify an identity value?

Do you agree that we need transform_* versions of algorithms?

Ranges plan papers (P2214R2 for C++23, P2760R1 for C++26) express a preference for views over algorithms.

This would call for ONLY reduce, inclusive_scan, and exclusive_scan, because users can get the functionality of their transform_* versions with views::transform and views::zip_transform.

We want transform_* versions of algorithms for performance, not functionality reasons. Does SG1 agree? This means:

• Prefer unary transform_reduce(range, 0.0, plus{}, get_element<0>{}) over
  • reduce(views::transform(get_element<0>{}, range), 0.0, plus{}) or
  • reduce(range, 0.0, [] (auto tup1, auto tup2) { return get<0>(tup2) + get<0>(tup2); })
• Prefer binary transform_reduce(r1, r2, T{}, plus{}, multiplies{}) over
  • reduce(views::zip_transform(multiplies{}, r1, r2), T{}, plus{})

Do you agree that we need a way for users to specify an identity value?

• Knowing the identity of a binary operator makes parallelizing reductions easier
• It depends on both a binary operator and the reduction result type
• Need not be the same as the initial value; interface must not confuse the two

Lack of an identity value complicates parallelization

Identity value != initial value

C++17 std::reduce takes an optional initial value T init.

std::vector<int> v{5, 7, 11};
const int init = 3;
auto result = std::reduce(v.begin(), v.end(), init, std::plus{});
assert(result == 26); // 3 + 5 + 7 + 11

Initial value defaults to typename std::iterator_traits::value_type{}, which just happens to be identity for plus on arithmetic types, though the algorithm doesn't need an identity.

Identity value id, if it exists, satisfies op(x, id) equals x and op(id, x) equals x for all x.

Identity need not be T{} for all operators and types

• For std::multiplies{} it's T(1)
• For "addition" in the max-plus ("tropical") algebra it's -Inf

Lack of identity hinders optimizing parallel reductions

Easy way to parallelize reduce:

1. Partition range into one contiguous subrange per processor
2. Each processor initializes its accumulator to identity (e.g., zero)
3. Use binary operator to combine accumulator with next element and store in accumulator
4. Combine all processors' accumulator values in some order

What if implementation has no identity value?

1. Same
2. Each processor initializes its accumulator to first element of its subrange
3. Same (starting with second element)
4. Same

This is functionally equivalent to calling reduce_first on each processor's subrange, instead of reduce. Why is this bad?

• Destroys overalignment of each subrange
  • SIMD-ization more complicated
  • May hinder other optimizations (e.g., some memcpy acceleration expects overalignment)
  • Users probably went through trouble to ensure overalignment
• Per-processor code gets longer and more branchy

Other parallel programming models let users specify identity

• Draft Fortran 2023 Standard (REDUCE clause)

• OpenMP ("used as the initializer for private copies of reduction list items")

• Kokkos (user-defined init(value_type& value))

• oneTBB (parallel_reduce parameter)

• SYCL (specialize sycl::known_identity class template for a custom reduction operation)

• NumPy: Can optionally specify identity of ufunc ("universal function," an elementwise binary operation); reductions use initial value or work like reduce_first if identity not specified

Use cases for default identity value

// Should this even work?
// Deduce binary op: std::plus{}
// Deduce identity: range_value_t<R>{}
auto result1 = std::ranges::reduce(exec_policy, range);

// Deduce identity: range_value_t<R>{}
auto result2 = std::ranges::reduce(exec_policy, range, std::plus{});

// If range_value_t<R> is arithmetic or std::complex,
// deduce identity: range_value_t<R>(1)
auto result3 = std::ranges::reduce(exec_policy, range, std::multiplies{});

// Should this even work?
// The "identity" is really -Inf if range_value_t<R> has that,
// but this is one of those cases where users might not like that behavior.
auto result4 = std::ranges::reduce(exec_policy, range, std::ranges::min);

Interface sketch for specifying identity

Design 1: reduce_identity<T>{value}

template<semiregular T>
struct reduce_identity {
  T value{};
};

• Default identity value is T{}
• Users have two ways to provide nondefault value
  1. reduce_identity{nondefault_value}
  2. Specialize reduce_identity<T> so declval<reduce_identity<T>>().value is the value

For example, users can inherit specialization from constant_wrapper.

namespace impl {
  inline constexpr my_number some_value = /* value goes here */;
}
template<class T>
struct reduce_identity<my_number> : 
  constant_wrapper<impl::some_value>
{};

Use cases:

// User explicitly opts into "most negative integer" as the identity for min.
// This should not be the default, as the C++ Standard Library has no way
// to know whether this represents a valid input value.
constexpr auto lowest = std::numeric_limits<int>::lowest();
auto result5 = std::ranges::reduce(exec_policy, range,
  std::ranges::min, reduce_identity{lowest});

// range_value_t<R> is float, but identity value is double
// (even though it's otherwise the default value, zero).
// std::plus<void> should use operator()(double, double) -> double
auto result6 = std::ranges::reduce(exec_policy, range,
  std::plus{}, reduce_identity{0.0});

Advantages:

• Algorithms can overload on it without risk of ambiguity
• Aggregate: minimizes compiler effort
• Does not impose extra requirements on binary function

Disadvantages:

• Algorithm can't use this to deduce default identity from binary operation
• User specialization of reduce_identity<T> will take effect for all binary operations on T

Design 2: reduce_operation{binary_op, value}

• Encapsulate binary operation and identity into a single argument
• Make it easier for implementation to deduce default identity value
• User specialization only affects that (operation, value type) pair
• It's still a struct, so algorithms can overload on it without ambiguity

auto result7 = std::ranges::reduce(exec_policy, range,
  reduce_operation{custom_binary_op, custom_value});

constexpr auto minus_Inf = -std::numeric_limits<float>::infinity();
auto result5 = std::ranges::reduce(exec_policy, range,
  reduce_operation{std::ranges::min, minus_Inf};

// range_value_t<R> is float
auto result6 = std::ranges::reduce(exec_policy, range,
  reduce_operation{std::plus{}, /* double */ 0.0};

// Deduce identity: range_value_t<R>{}
auto result2 = std::ranges::reduce(exec_policy, range, std::plus{});

Here is an implementation sketch.

// Keep this an aggregate if possible.
// We can constrain BinaryOp on being invocable with T, T,
// but we don't know range_value_t<R> at this point,
// so we don't know yet if the algorithm is well-formed with BinaryOp.
template<class BinaryOp, class T = invoke_result<BinaryOp, T, T>>
  requires semiregular<T> && invocable<BinaryOp, T, T>
struct reduce_operation {
  BinaryOp op;
  T value{};
};

// std::plus with arithmetic types -> T{}.
//
// We can't specialize for all T, because users might have
// defined std::plus<T> such that T{} is not the identity.
template<class U, class T>
  requires(
    (is_arithmetic_v<T>) &&
    (is_void_v<U> || (is_arithmetic_v<U> && is_convertible_v<U, T>))
  )
struct reduce_operation<std::plus<U>, T> {
  static constexpr std::plus<U> op{};
  static constexpr T value = T{};
};

// (std::plus or std::multiplies) with std::complex takes more effort,
// since e.g., std::complex<double> + float is not well-formed.

// std::multiplies with arithmetic types -> T(1).
//
// We can't specialize for all T, because users might have
// defined std::multiplies<T> such that T(1) is not the identity.
template<class T, class U>
  requires(
    (is_arithmetic_v<T>) &&
    (is_void_v<U> || (is_arithmetic_v<U> && is_convertible_v<U, T>))
  )
struct reduce_operation<T, std::multiplies<U>> {
  static constexpr std::multiplies<U> op{};
  static constexpr T value = T(1);
};

// ... more specializations for Standard binary function objects ...

Some design issues:

1. Users can't override default identity values for known binary operators and reduction result types.
   • (Feature, not a bug?)
2. BinaryOp must be copy-constructible or move-constructible.
   • That's fine for parallel algorithms, but not strictly necessary for non-parallel algorithms.
   • On the other hand, if you want non-copyable BinaryOp, maybe you should use fold_* instead.

*transform_view not trivially copyable even if function is

Implementations on accelerators with their own memory (e.g., GPU) want to memcpy (or equivalent) the algorithm's arguments to accelerator memory. This only works if the arguments are all trivially copyable.

Problem: transform_view and zip_transform_view might not necessarily be trivially copyable, even if their function object is.

• Example: lambda with by-copy capture of an int
• Probably a wording bug (relates to movable-box)

This may prevent use of accelerators for such views.

(Ranges) views generally can hinder optimizations

Accelerator vendors have reasons not to reimplement views

Accelerator-based implementations generally do not want to reimplement the whole Standard Library. They want to work with any portable C++ Standard library, and with its implementations of views like *transform_view.

• P2500 expresses our end goal: accelerator vendors can customize algorithms in a Standard Library implementation that they didn't write.

  • This expresses how accelerators are used: as targeted optimizations for specific functions, not for the whole application.

• Users may compile mostly with host compiler (e.g., Clang, GCC, MSVC) and use a special compiler for accelerators.

  • Accelerator-based implementations then reimplement (and/or specialize for custom execution policies) only a few algorithms, not views.

  • Accelerator vendors don't want to force users to adopt a new Standard Library implementation: it would break ABI and prevent them from using host-compiler-specific extensions.

• Users may assemble the same range to run on either host or accelerator.

Ranges are not designed for algorithm customization

• Views like *transform_view do not generally give access to their members.

  • Ranges designers might consider this a feature, not a bug.

  • Many ranges are move-only; *transform_view holds them.

  • Contrast with std::execution, where sender customization is the whole point.

• As a result, parallel reduce can't specialize on *transform_view range input in ways that need to take apart the view.

  • e.g., memcpy just the function and range, not the view, to the accelerator's memory

  • e.g., translate reduce(par, views::zip_transform(op1, r1, r2), op2) to a binary transform_reduce(par, r1, r2, op1, op2) internally

• This affects other views too; e.g., can't transform ranges::for_each of cartesian_product_view of iota_views into a custom library function call for multidimensional loops.

Design implications

• Provide transform_* versions of algorithms, rather than relying on transform_view
• Provide binary transform_reduce, instead of relying on zip_transform_view

Expose built-in arithmetic on arithmetic types, where possible

Existing parallel programming models prefer this

• Vendors like to use existing parallel programming models to implement C++ parallel algorithms
  • Someone has already optimized them
  • Some integrate with the compiler (e.g., OpenACC, OpenMP) and thus have more knowledge about code than a library could
• Some parallel programming models do not support custom reduction operations
• Built-in types and arithmetic can exploit more reduce-accelerating hardware features
  • atomic +=
  • SIMD reductions
  • Network hardware for reductions

Some design implications

• Prefer unary transform_reduce(range, T{}, 0.0, plus{}, get_element<0>{}) over
  • reduce(views::transform(get_element<0>{}, range), 0.0, plus{}) or
  • reduce(range, 0.0, [] (auto tup1, auto tup2) { return get<0>(tup2) + get<0>(tup2); })
• Prefer binary transform_reduce(r1, r2, T{}, plus{}, multiplies{}) over
  • reduce(views::zip_transform(multiplies{}, r1, r2), T{}, plus{})
• Consider projections, perhaps even with binary transform_reduce, where they make the interface hard to use (four function parameters in a row)