# Reductions: current implementation and known problems/TODOs

**NOTE**: This document is a quick draft. It is written to help developers of SYCL headers/library to understand the current status, currently used algorithms and known problems.

---
---
# Implementation details: `reduction` in `parallel_for()` accepting `nd_range`

## Reduction inside 1 work-group - the main building block for reduction

The reduction algorithm for 1 work-group depends currently on combination/availability of 2 features:
- `fast atomics` (i.e. operations available on the target device), such as fetch_add(), fetch_min(), etc.
- `fast reduce` operation for work-group, i.e. `reduce_over_group()`

So, if the reduction operation/type has both `fast atomics` and `fast reduce`, then the reduction on work-group with 1 reduction variable does the following: the elements inside the work-group are reduced using `ext::oneapi::reduce` and the final result is atomically added to the final/global reduction variable.

```c++
    // Call user's lambda function or functor. Reducer.MValue gets initialized there.
    typename Reduction::reducer_type Reducer;
    KernelFunc(NDIt, Reducer);

    typename Reduction::binary_operation BOp;
    // Compute the partial sum for the work-group.
    // Then the 0-th work-iteam stores it to the final/global reduction variable.
    Reducer.MValue = ext::oneapi::reduce(NDIt.get_group(), Reducer.MValue, BOp);
    if (NDIt.get_local_linear_id() == 0)
      Reducer.atomic_combine(Reduction::getOutPointer(Out));
```

The most general case is when the `fast atomics` and `fast reduce` are not available.
It computes the partial sum using tree-reduction loop and stores the partial sum for the work-group to a global array of partial sums, which later must be also reduced (by using additional kernel(s)). This general case algorithm also requires allocation of a local memory for the tree-reduction loop:
```c++
    // Call user's functions. Reducer.MValue gets initialized there.
    typename Reduction::reducer_type Reducer(ReduIdentity, BOp);
    KernelFunc(NDIt, Reducer);

    size_t WGSize = NDIt.get_local_range().size();
    size_t LID = NDIt.get_local_linear_id();
    // Copy the element to local memory to prepare it for tree-reduction.
    LocalReds[LID] = Reducer.MValue;
    if (!IsPow2WG)
      LocalReds[WGSize] = ReduIdentity;
    NDIt.barrier();

    // Tree-reduction: reduce the local array LocalReds[:] to LocalReds[0]
    // LocalReds[WGSize] accumulates last/odd elements when the step
    // of tree-reduction loop is not even.
    size_t PrevStep = WGSize;
    for (size_t CurStep = PrevStep >> 1; CurStep > 0; CurStep >>= 1) {
      if (LID < CurStep)
        LocalReds[LID] = BOp(LocalReds[LID], LocalReds[LID + CurStep]);
      else if (!IsPow2WG && LID == CurStep && (PrevStep & 0x1))
        LocalReds[WGSize] = BOp(LocalReds[WGSize], LocalReds[PrevStep - 1]);
      NDIt.barrier();
      PrevStep = CurStep;
    }

    // Compute the partial sum/reduction for the work-group and store it to the global array of partial sums.
    if (LID == 0) {
      size_t GrID = NDIt.get_group_linear_id();
      typename Reduction::result_type PSum =
          IsPow2WG ? LocalReds[0] : BOp(LocalReds[0], LocalReds[WGSize]);
      if (IsUpdateOfUserVar) // depends on accessor read/write mode or on property::reduction::initialize_to_identity
        PSum = BOp(*(Reduction::getOutPointer(Out)), PSum);
      Reduction::getOutPointer(Out)[GrID] = PSum;
    }
```

The *additional kernel(s)* that reduce the partial sums look very similar to the *main* kernel, the only difference is that they do not call user's function.

---
TODO #1 (Performance): After the reduction accepting 'range' was implemented it became possible to simply re-use that implementation for the *additional kernel* runs. Instead of running the *additional kernel* several times (if the number of partial sums is still big), it is better to run the parallel_for accepting 'range' once.

---
TODO #2 (Performance): There are 4 implementations for `parallel_for()` accepting nd_range and 1 reduction variable now:
*  `fast atomics` and `fast reduce` (A)
*  `fast atomics` and `no fast reduce` (B)
*  `no fast atomics` and `fast reduce` (C)
*  `no fast atomics` and `no fast reduce` (D)

Accordingly to the latest performance results/experiments it may be more efficient to replace (A) and (B) with one implementation using `fast atomics` for 1) reducing inside work group - use atomic ops to 1 scalar allocated in group local memory, 2) use atomics to store the partial sum for work-group to final/global reduction variable.

---
TODO #3 (Performance): Currently, there is more or less efficient implementation for 1 reduction variable. Using `fast atomics` and `fast reduce` does give better results. The implementation for any number of reduction variables used in same `parallel_for()` currently uses only the most basic implementation using `tree-reduction` algorithm for all reduction operations/types even those that do have `fast atomics` or `fast reduce` algorithm.

---
---
# Implementation details: `reduction` in `parallel_for()` accepting `range`

For `parallel_for()` accepting 1 reduction variable the implementation chooses one of 3 algorithms depending on existence of `fast atomics` and `fast reduce` features in the target device, which gives 3 variants:
*  `fast atomics` (A)
*  `no fast atomics` and `fast reduce` (B)
*  `no fast atomics` and `no fast reduce` (C)

The implementation of reduction accepting `range` has more freedom comparing to `nd_range` case as the order of reducing elements in each work-items is not specificed and the work-group sizes not specified too.

The implementation queries the target device on the number of execution units and max work-group-size, after which it chooses number and size of work-groups, then calls `parallel_for()` accepting `nd_range` inside the original/user's `parallel_for()` accepting `range`. Each of work-items in that nd-range space gets 1 or more indices from the global index space for which it calls user's function.

For the case (A) `fast atomics` the implementation looks this way:

```c++
  // Note that this parallel_for() is not user's parallel_for(). It is called inside user's one.
  CGH.parallel_for<Name>(NDRange, [=](nd_item<1> NDId) {
    // Call user's functions. Reducer.MValue gets initialized there.
    typename Reduction::reducer_type Reducer;
    reductionLoop(Range, Reducer, NDId, KernelFunc); // Each work-item handles 1 or many indices
                                                     // from the original/user's global index space.

    // Store the accumulated partial sum for the work-item to local var holding the partial sum for the work-group.
    auto LID = NDId.get_local_id(0);
    if (LID == 0)
      GroupSum[0] = Reducer.getIdentity();
    sycl::detail::workGroupBarrier();
    Reducer.template atomic_combine<access::address_space::local_space>(
        &GroupSum[0]);

    // Store the partial sum for the work-group to the final/global reduction variable 
    sycl::detail::workGroupBarrier();
    if (LID == 0) {
      Reducer.MValue = GroupSum[0];
      Reducer.template atomic_combine(Reduction::getOutPointer(Out));
    }
  });
```
Variants (B) and (C) use the same approach. The only difference is how the partial sums accumulated per work-items are combined into a partial sum for the work-group. Instead of fast atomic to a local scalar either the `fast reduce` or `tree-reduce-loop` used there. 

---

TODO #4 (Performance): Some devices may provide unique-thread-id where the number of worker threads running simultaneously is limited. Such feature opens way for more efficient implementations (up to 2x faster, especially on many stacks/tiles devices). See this extension for reference: https://github.com/intel/llvm/pull/4747

---
---

# Missing `reduction` functional features to be implemented soon:

### 1) Support `sycl::queue::parallel_for()` shortcuts accepting reduction variables.
Currently, only `sycl::queue::parallel_for()` accepting `nd_range` and 1 reduction variable is supported.
The rest of this work is temporarily blocked by XPTI instrumentation that need the last argument to be hidden from user and having a default value, which is not possible if `parallel_for` starts accepting parameter packs like this one:
```c++
 // Parameter pack acts as-if: Reductions&&... reductions, const KernelType &kernelFunc
 template <typename KernelName, int Dims, typename... Rest>
 event parallel_for(range<Dims> numWorkItems, Rest&&... rest);
```

The problem is known, the fix in SYCL headers is implemented: https://github.com/intel/llvm/pull/4352 and is waiting for some re-work in XPTI component that must be done before the fix merge.

---
### 2) Support identity-less reductions even when the reduction cannot be determinted automatically.

Currently identity-less reductions are supported, but only in cases when sycl::has_known_identity<BinaryOperation, ElementType> returns true.
When sycl::has_known_identity returns false, the implementation of the reduction may be less efficient, but still be functional.

---
## Note the document is not finished. It should be updated with more details soon.