---
title: "Crushing CPUs with 879 GB/s Reductions in CUDA"
date: 2022-01-28
description: "GPU code beats optimized CPU parallel reductions by 10x, reaching 879 GB/s. CUB achieves 94% bandwidth saturation while CPU barely hits 60%."
tags: [Less Slow]
source: https://ashvardanian.com/posts/cuda-parallel-reductions/
author: Ash Vardanian
---



GPU acceleration [can be trivial](https://unum.cloud/blog/2022-01-26-cupy/) for Python users.
Follow CUDA installation steps carefully, replace `import numpy as np` with `import cupy as np`, and you will often get the 100x performance boosts without breaking a sweat.
Every time you write magical one-liners, remember a systems engineer is making your dreams come true.

---

A couple of years ago, when I was giving a [talk](https://www.youtube.com/watch?v=AA4RI6o0h1U) on the breadth of GPGPU technologies, I published a repo.
A repo with various cool GPGPU staff in [Vulkan](https://www.vulkan.org), [Halide](https://halide-lang.org), and most importantly, __8x__ implementations of `std::accumulate` in OpenCL.
I know what you are thinking:

> Eight ways to sum numbers?! Are you 🌰s?!

Don't worry!
We are now back with more!
No more OpenCL, no more [SyCL](https://www.khronos.org/sycl/).
__This time we focus on [Parallel STL](#parallel-algorithms), [SIMD](#simd-avx2), [OpenMP](#openmp) and, of course, [CUDA](#cuda-with-shfl), [CUB](#cub) and [Thrust](#thrust)!__

> You can find all the sources [hosted on our __GitHub__](https://github.com/ashvardanian/ParallelReductions).
> Benchmarks were run with the newest stable versions of software available at the time: Ubuntu 20.04, GCC 10.3, CUDA Toolkit 11.6, Thrust 1.15, oneTBB 2021.5 and TaskFlow 3.3.

---

> This article went through 2 updates.
> First, the [OpenMP numbers](#update-on-openmp) were corrected.
> Later, I reflected on the results and came to a shocking conclusion.
> It thematically belongs in the middle, but I [published it separately](/posts/ddr4-bandwidth/) to avoid spoilers and keep it chronological.

## C++ and STL

The canonical serial solution for this problem in C and C++ would be:

```cpp
float sum = 0
for (; begin != end; ++begin)
    sum += *begin;
```

For simple tasks like this, there is also an STL version.
Let's put a sample an example with 1 GB worth of `float` numbers:

```cpp
std::vector<float> numbers(1024 * 1024 * 256);
std::fill(numbers.begin(), numbers.end(), 1.f);
auto sum = std::accumulate(numbers.begin(), numbers.end(), 0.f);
```

Time for some trivia questions. What will be the sum?
Due to rounding errors, the result of sequential accumulation will be significantly less accurate than doing it in batches or parallel tree-like reductions.
To avoid it, you would need more memory.
One more `float` for the compensation part, as in the [Kahan summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) algorithm.
Or simply using a `double` for the accumulation.
The latter being equally fast on modern x86 chips and easier to implement:

- __5.2 GB/s__ when accumulating into `float`, with a 93% error.
- __5.3 GB/s__ when accumulating into `double`, with a 0% error.

## SIMD: AVX2

We also implemented three [AVX2 SIMD](https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#techs=AVX,AVX2) variants for intra-thread acceleration:

1. Naive 8x lane `float` accumulation.
2. Kahans 8x lane `float` accumulation.
3. Conversions and 4x lane `double` accumulation.

They all utilize heavy instructions, so some downclocking occurs, but performance is excellent, and rounding errors disappear in the last two variants.
The (1) and (3) SIMD code is trivial.
We will only post the Kahan version for clarity.

```cpp
float operator()(float const * begin, float const * end) noexcept {
    auto it = begin;

    // SIMD-parallel summation stage
    auto sums = _mm256_set1_ps(0);
    auto compensations = _mm256_set1_ps(0);
    auto t = _mm256_set1_ps(0);
    auto y = _mm256_set1_ps(0);
    for (; it + 8 < end; it += 8) {
        y = _mm256_sub_ps(_mm256_loadu_ps(it), compensations);
        t = _mm256_add_ps(sums, y);
        compensations = _mm256_sub_ps(_mm256_sub_ps(t, sums), y);
        sums = t;
    }

    // Cross-lane horizontal reduction
    sums = _mm256_hadd_ps(sums, _mm256_permute2f128_ps(sums, sums, 1));
    sums = _mm256_hadd_ps(sums, sums);
    sums = _mm256_hadd_ps(sums, sums);
    auto sum = _mm256_cvtss_f32(sums);

    // Serial summation of the remainder
    for (; it != end; ++it)
        sum += *it;

    return sum;
}
```

Results on a single CPU core:

1. __22.3 GB/s__ naively accumulating `float`, with a 50% error.
2. __10.9 GB/s__ accumulation with Kahans method, with 0% error.
3. __17.0 GB/s__ accumulating as `double`s, with 0% error.

## OpenMP

Not to be confused with OpenMPI, [Open Multi-Processing](https://en.wikipedia.org/wiki/OpenMP) is probably the oldest Multi-Threading (not Multi-Processing!) standard.
It has run on anything from desktops to supercomputers since 1997 and is widely supported by Fortran, C and C++ compilers.
Aside from the oldest `#pragma`s, they also support parallel reductions:

```cpp
float operator()(float const * begin, float const * end) noexcept {
    float sum = 0;
    size_t const n = end - begin;
#pragma omp parallel for default(shared) reduction(+ : sum)
    for (size_t i = 0; i != n; i++)
        sum += begin[i];
    return sum;
}
```

We tried [every OpenMP reduction tutorial](https://www.sciencedirect.com/topics/computer-science/parallel-reduction) but couldn't make it work faster than the most basic serial version of `std::accumulate`.
Result: __5.4 GB/s__.

### Update on OpenMP

After a few recommendations on [Reddit](https://www.reddit.com/r/cpp/comments/se2eks/best_cpu_vs_gpu_879_gbs_reductions_in_c/?utm_source=share&utm_medium=web2x&context=3), an issue with compilation flags was resolved.
At 100% CPU utilization OpenMP scored __51.5 GB/s__.
After that, we disabled dynamic scheduling with `omp_set_dynamic(0)` and reduced the number of threads, reaching the result of __80 GB/s__.

## Parallel Algorithms

STL ships with [`std::execution`](https://en.cppreference.com/w/cpp/algorithm/execution_policy_tag_t) policies since 17th edition.
Hope being, that multi-threaded `<algorithm>`s are easy to use and at least as good as average parallel code.

```cpp
auto sum = std::reduce(std::execution::par_unseq, numbers.begin(), numbers.end(), 0.f);
```

Result: __5.3 GB/s__.
Wait, it's the same as we got for single-threaded code.
We forgot that GCC relies on Intel's [Thread Building Blocks](https://www.intel.com/content/www/us/en/develop/documentation/oneapi-programming-guide/top/api-based-programming/intel-oneapi-threading-building-blocks-onetbb.html) to implement parallel algorithms.
Let's update our `CMakeLists.txt`:

```cmake
FetchContent_Declare(
    TBB
    GIT_REPOSITORY https://github.com/oneapi-src/oneTBB.git
    GIT_TAG        v2021.5.0
)
FetchContent_Populate(TBB)
include_directories(BEFORE ${TBB_SOURCE_DIR})
target_link_libraries(reduce_bench TBB::tbb)
```

Run again and hurray!

- __80 GB/s__ in `std::execution::par` reductions.
- __87 GB/s__ in `std::execution::par_unseq` reductions.

Now we are going somewhere!
Can we go there faster?

## SIMD + Threads

If we take our AVX2 SIMD (#3) implementation and spawn 64x `std::threads`, whenever a new task comes, we can still go a little faster: __89 GB/s__.
In the best-case scenario, if we always had a thread-pool around, with ~20 idle threads, we could do even better, up to 200 GB/s.

![Fastest CPU-RAM Links](/cuda-parallel-reductions/ram_bandwidths.png)

Not more.
Not until [we switch](/posts/server-supercycle) from DDR4 to DDR5 and from 8-channel to 12-channel RAM.
Extrapolating the [Hash-Table benchmarks](/posts/apple-m1) of Apple M1 Max laptops, we can expect massive boosts in DDR5-powered servers.
Sapphire Rapids may even get in-package High Bandwidth Memory, but for now, only GPUs have that!

## Switching to GPUs

It would have been interesting to run the same OpenCL scripts on both CPUs and GPUs.
Unfortunately, we lost that opportunity about two years ago, when AMD [dropped support for OpenCL](https://www.reddit.com/r/Amd/comments/aun3y2/amd_cpus_no_longer_have_opencl_support/) in their CPU lineup.
As we are only limited to a GPU, we will skip OpenCL this time, assuming that, on average, its kernels are at least 30% slower than CUDA.

| GPU Model                                                                  | Size  |  Type  |   Bus    | Bandwidth  |
| :------------------------------------------------------------------------- | :---: | :----: | :------: | :--------: |
| [3090](https://www.techpowerup.com/gpu-specs/geforce-rtx-3090.c3622)       | 24 GB | GDDR6X | 384 bit  |  936 GB/s  |
| [3090 Ti](https://www.techpowerup.com/gpu-specs/geforce-rtx-3090-ti.c3829) | 24 GB | GDDR6X | 384 bit  | 1'018 GB/s |
| [A100](https://www.techpowerup.com/gpu-specs/a100-sxm4-80-gb.c3746) SXM4   | 80 GB | HBM2e  | 5120 bit | 2'039 GB/s |

Before starting our experiment, let's determine our upper bound.
The A100 is a datacenter GPU, and 3090 Ti hasn't reached the market yet, so we won't get 1 TB/s this time around.

## CUDA with SHFL

There is the old way "to CUDA" and the new way.
The old one uses just `__shared__` memory.
The post-Kepler method recommends using shuffles for thread scheduling.

```cpp
__inline__ __device__ float cu_reduce_warp(float val) {
    val += __shfl_down_sync(0xffffffff, val, 16);
    val += __shfl_down_sync(0xffffffff, val, 8);
    val += __shfl_down_sync(0xffffffff, val, 4);
    val += __shfl_down_sync(0xffffffff, val, 2);
    val += __shfl_down_sync(0xffffffff, val, 1);
    return val;
}

__global__ void cu_reduce_warps(float const *inputs, unsigned int input_size, float *outputs) {
    float sum = 0;
    for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
            i < input_size;
            i += blockDim.x * gridDim.x)
        sum += inputs[i];

    __shared__ float shared[32];
    unsigned int lane = threadIdx.x % warpSize;
    unsigned int wid = threadIdx.x / warpSize;

    sum = cu_reduce_warp(sum);
    if (lane == 0)
        shared[wid] = sum;

    // Wait for all partial reductions
    __syncthreads();

    sum = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
    if (wid == 0)
        sum = cu_reduce_warp(sum);

    if (threadIdx.x == 0)
        outputs[blockIdx.x] = sum;
}
```

This sample more or less repeats the CUDA tutorial but replaces the [deprecated](https://developer.nvidia.com/blog/using-cuda-warp-level-primitives/) `__shfl_down` with `__shfl_down_sync` in the first function.
Mouthful, compared to `std::accumulate`, but it gets the work done: __817 GB/s__ 💣💣

## Thrust

Thrust might just be my favorite high-level library.
Generic enough to be considered an STL extension and full of intrguing technical solutions.
If you don't want to write a custom kernel for a GPU, you can probably compose a multi-round alternative with Thrust.

```cpp
thrust::reduce(numbers.begin(), numbers.end(), float(0), thrust::plus<float>());
```

The result: __743 GB/s__.
A 9% reduction, compared to the hand-written kernel, but hardly slow.

Those of us waiting for NVCC implementation of Parallel Algorithms will use it indirectly.
Others can use it directly or go one ~~paragraph~~ layer lower to CUB.

## CUB

CUBs main benefit is stricter control over memory allocations.
Parallel and concurrent algorithms often require extra memory compared to serial analogues.
To keep the interface familiar and straightforward, Thrust may allocate temporary memory internally.

With CUB, you manage memory yourself and control the algorithm selection more explicitly, resulting in faster code but a much heavier codebase.
For example, you would typically call the same function twice.
First time only to estimate the needed temporary memory capacity.

```cpp
float operator()(float const *b, float const *e) {

    // Reuse the following arrays
    thrust::device_vector<float> gpu_inputs(b, e);
    thrust::device_vector<float> gpu_sums(1);
    thrust::host_vector<float> cpu_sums(1);
    thrust::device_vector<uint8_t> temporary;

    // CUB can't handle large arrays with over 2 billion elements!
    assert(gpu_inputs.size() < std::numeric_limits<int>::max());
    auto num_items = static_cast<int>(gpu_inputs.size());
    auto ins = gpu_inputs.data().get();
    auto outs = gpu_sums.data().get();

    // Determine temporary device storage requirements
    cudaError_t error;
    void *temps = nullptr;
    size_t temps_size = 0;
    error = cub::DeviceReduce::Sum(temps, temps_size, ins, outs, num_items);
    assert(error == cudaSuccess);

    // Allocate temporary storage, if needed
    if (temps_size > temporary.size())
        temporary.resize(temps_size);
    temps = temporary.data().get();

    error = cub::DeviceReduce::Sum(temps, temps_size, ins, outs, num_items);
    assert(error == cudaSuccess);
    cudaDeviceSynchronize();
    cpu_sums = gpu_sums;
    return cpu_sums[0];
}
```

CUB provides shortcuts for the most common binary commutative operators: sum, min, max, argmin, argmax.
So instead of invoking [`cub::DeviceReduce::Reduce`](https://nvlabs.github.io/cub/structcub_1_1_device_reduce.html) in this case, we run `cub::DeviceReduce::Sum`.

Result: __879 GB/s__ 🔥🔥🔥

CUB is an excellent library by itself - a huge productivity booster.
Especially the device-scale functions, like `DeviceHistogram`, `DevicePartition`, `DeviceRadixSort`, `DeviceRunLengthEncode`, `DeviceScan`.
Our internal on-GPU compression library, for example, has a multi-stage pipeline, where parts were initially written with Thrust.
Then accelerated with CUB and only in most important places, rewritten in raw CUDA.

## Tensor Cores

This is one of those cases where [Tensor Cores](https://developer.nvidia.com/blog/programming-tensor-cores-cuda-9/) can't help us much.
We can program the `wmma::` intrinsics in CUDA directly.
In that case, we can reshape the input array into a very long matrix and multiply it by small auto-generated matrices on the fly.
Every `wmma::fragment<wmma::matrix_b, 16, 16, 8, wmma::precision::tf32, wmma::row_major>` would contain mostly zeros, except the first column of ones.
Every warp would perform its reduction, accumulating 16²=256 numbers into a column of 16 sums.
One per row.
It's a lot of wasted compute capacity, but given the insane efficiency of Tensor Cores, we could have gotten a slight improvement.

Still, assuming we have already reached __94% of theoretical memory bandwidth__, we are definitely in the red zone of diminishing returns.

## Roundup

We repeated the benchmarks for every power-of-two array size between 32 MB and 4 GB.
If an image is worth a thousand words to you, enjoy!

![Parallel Reduction Speeds](/cuda-parallel-reductions/plot.png)

Again, if you are curious to see how your CPU/GPU would compete, here is the [GitHub](https://github.com/ashvardanian/ParallelReductions) link for these benchmarks.
It should be pretty easy to install if you have ever used CMake.

Intel is preparing its [Ponte Vecchio](/posts/server-supercycle/#intel) GPU launch.
And AMD is trying to attract more customers for their [MI200](/posts/server-supercycle/#amd) GPUs.
With all that action, we should expect more first-party HPC libraries, like Thrust and CUB.
Before then, we will keep reinventing the wheel, often implementing stuff from scratch and writing about it occasionally.
All to accelerate AI research and build the [__fastest persistent DBMS__](https://unum.cloud/blog/2021-11-25-ycsb/) the world as ever seen 😉