As my first foray into CUDA programming I set myself a challenge to experiment with different methods of finding the total sum of a large array of numbers, trying different strategies exposed via CUDA to see how different approaches compare and what comes out fastest. A lot of my GPU experience previously on PC has been through DirectX 12 or Vulkan so this felt like a good opportunity to see how programming CUDA compares for a non-graphics or games related workload.

Source Code: https://github.com/comeradealexi/CUDAParallelAddBenchmarks/blob/main/kernel.cu

Test Setup

Each test reads from a large array of integers, adding every entry together and returning the sum. Given the parallel nature of GPUs, we use the atomicAdd function to allow separate threads and blocks to increment the final output value.

For the various test scenarios which are detailed below, we will run with a suite of different block sizes as well as the number of values each thread processes from the input array in an attempt to find the most optimal setup. Block sizes will differ from 32 up to 1024 which has an effect on occupancy as well as the tests involving shared memory.

My expectation is that we’ll get faster results as we reduce the number of times we call atomicAdd, however this isn’t quite what my results find.

With the specific test scenarios below, we will run with the following setups:

• 512MiB Buffer & 1 Kernel Run
• 4 MiB Buffer & 64/128 Kernel Runs

I wanted to test the difference in the results between a very large buffer as well as a small kernel accessed multiple times to ensure I observe the difference in the memory speed overhead. I’m expecting that the 4MiB buffer tests, ran multiple times will remain in the GPU caches and better tests how the code is performing, while I expect the very large buffers to have memory access overheads, potentially resulting in differing results.

Note: In order to ensure stable GPU frequencies, these tests have been run with the following enabled: nvidia-smi -lgc 2800

Test Scenarios

template<unsigned int VALUES_PER_THREAD>
__global__ void Add_AtomicPerValuePerThread(int* output, const int* input)
{
    const unsigned int workIndex = (threadIdx.x + (blockDim.x * blockIdx.x)) * VALUES_PER_THREAD;

    for (unsigned int i = 0; i < VALUES_PER_THREAD; i++)
    {
        int value = input[workIndex + i];
        atomicAdd(&output[0], value);
    }
}

template<unsigned int VALUES_PER_THREAD>
__global__ void Add_AtomicPerThread(int* output, const int* input)
{
    const unsigned int workIndex = (threadIdx.x + (blockDim.x * blockIdx.x)) * VALUES_PER_THREAD;

    int total = 0;    
    for (unsigned int i = 0; i < VALUES_PER_THREAD; i++)
    {
        total += input[workIndex + i];
    }
    atomicAdd(&output[0], total);
}

template<unsigned int VALUES_PER_THREAD>
__global__ void Add_ReduceAddSyncPerWarp(int* output, const int* input)
{
   const unsigned int workIndex = (threadIdx.x + (blockDim.x * blockIdx.x)) * VALUES_PER_THREAD;

    int total = 0;
    for (unsigned int i = 0; i < VALUES_PER_THREAD; i++)
    {
        total += input[workIndex + i];
    }

    // Sum of all 32 threads in a given warp
    int blockTotal = __reduce_add_sync(0xFFFFFFFF, total);

    // First of every 32 threads atomically increments final output
    if ((threadIdx.x & 31) == 0)
    {
        atomicAdd(&output[0], blockTotal);
    }
}

template<unsigned int VALUES_PER_THREAD>
__global__ void Add_AtomicSharedPerThread(int* output, const int* input)
{
    __shared__ int sharedData;
    if (threadIdx.x == 0)
    {
        sharedData = 0;
    }
    __syncthreads();

    const unsigned int workIndex = (threadIdx.x + (blockDim.x * blockIdx.x)) * VALUES_PER_THREAD;

    int total = 0;
    for (unsigned int i = 0; i < VALUES_PER_THREAD; i++)
    {
        total += input[workIndex + i];
    }

    atomicAdd(&sharedData, total);

    __syncthreads();

    if (threadIdx.x == 0)
    {
        atomicAdd(&output[0], sharedData);
    }
}

template<unsigned int VALUES_PER_THREAD>
__global__ void Add_SharedUniqueEntryPerThread(int* output, const int* input) // Global atomic increments = Total Threads / Block Size
{
    constexpr int SHARED_DATA_LENGTH = 1024;
    __shared__ int sharedData[SHARED_DATA_LENGTH];

    const unsigned int workIndex = (threadIdx.x + (blockDim.x * blockIdx.x)) * VALUES_PER_THREAD;

    int total = 0;
    for (unsigned int i = 0; i < VALUES_PER_THREAD; i++)
    {
        total += input[workIndex + i];
    }
    sharedData[threadIdx.x] = total;
    
    __syncthreads();

    if (threadIdx.x == 0) // First thread in the block loops over all shared data.
    {
        int blockTotal = 0;
        for (int i = 0; i < blockDim.x; i++)
        {
            blockTotal += sharedData[i];
        }
        atomicAdd(&output[0], blockTotal);
    }
}

template<unsigned int VALUES_PER_THREAD, unsigned int BLOCK_SIZE>
__global__ void Add_SharedWithReduceAddSyncPerWarp(int* output, const int* input) // Global atomic increments = Total Threads / Block Size
{
    constexpr int SHARED_DATA_LENGTH = (BLOCK_SIZE / 32) > 0 ? (BLOCK_SIZE / 32) : 1;
    __shared__ int sharedData[SHARED_DATA_LENGTH];

    const unsigned int workIndex = (threadIdx.x + (blockDim.x * blockIdx.x)) * VALUES_PER_THREAD;

    int total = 0;
    for (unsigned int i = 0; i < VALUES_PER_THREAD; i++)
    {
        total += input[workIndex + i];
    }

    int warpTotal = __reduce_add_sync(0xFFFFFFFF, total);

    if ((threadIdx.x & 31) == 0) // Each warp sums all threads in that warp.
    {
        // Divide by warp count
        sharedData[threadIdx.x / 32] = warpTotal;
    }
   
    __syncthreads();

    if (threadIdx.x == 0)
    {
        int blockTotal = 0;
        for (unsigned int i = 0; i < SHARED_DATA_LENGTH; i++)
        {
            blockTotal += sharedData[i];
        }

        atomicAdd(&output[0], blockTotal);
    }
}

Results

As you would expect, given the number of test scenarios, the differing combinations of block/grid sizes as well as our differing memory sizes, I’ve ended up with a lot of data. I’ll try to share it and analyse it in parts where the results really stood out to me.

Results have been collected on an Nvidia RTX 5060 running on Windows 11 with CUDA 13.2.1. Below are the device capabilities for the Nvidia RTX 5060:

Device 0:
  Name:                        NVIDIA GeForce RTX 5060
  WarpSize:                    32
  MaxThreadsPerBlock:          1024
  MaxGridSize(X):              2147483647
  MaxGridSize(Y):              65535
  MaxGridSize(Z):              65535
  MaxThreadsDim(X):            1024
  MaxThreadsDim(Y):            1024
  MaxThreadsDim(Z):            64
  SharedMemPerBlock:           49152
  SharedMemPerMultiprocessor:  102400
  MaxThreadsPerMultiProcessor: 1536
  MaxBlocksPerMultiProcessor:  24
  MultiProcessorCount:         30
  SMCount:                     30
  MinSMPartitionSize:          8
  SMCoscheduledAlignment:      8
  WQConcurrencyLimit:          8

512MiB x 1 Results:

Comparing tests with same block sizes:

Here we are comparing all 6 of our test scenarios with the same block and grid sizes, so each test issues exactly the same number of threads. All of these test scenarios achieve a potential occupancy of 100%.

The part that surprised me here is just how close majority of these results are in terms of performance, especially the first three results. Add_AtomicPerThread does 512 global atomic increments per block and it’s performing better than several tests which do fewer global atomic increments.Add_AtomicPerValuePerThread is the outlier here, with such a large number of atomic increments per thread, showing an 80% increase in performance compared to our fastest scenario.

How occupancy affects results:

Focusing in on the Add_AtomicPerThread scenario, it’s quite clear here how not achieving full occupancy is having a detrimental effect on the test results, to the point where the time is more than doubled when comparing 100% to 50% occupancy. This suggests that each SM is spending time idle, likely waiting for memory operations to complete and cannot spend that time processing more warps due to the tests not achieving full occupancy.

4MiB x 64 Results:

With the 4MiB tests, where each scenario is run 64 times back to back, we are expecting very different results to the single 512MiB tests because we expect the input memory to always remain within the GPU caches, so less time is spent waiting on memory operations. We’ll start by comparing the fastest results for each test scenario:

The results here align more with our initial expectations where the performance scales with the number of atomic increments per thread. But the difference is only small, comparing the fastest Add_AtomicSharedPerThread to to Add_AtomicPerThread, there are 512x more global atomic increments per thread, but performance is only 1.25x slower.

Importance of keeping the GPU saturated

The results below show the how increasing the amount of work done per-thread decreases performance. Each SM on the GPU is designed to have many blocks assigned to it, ready to execute when other threads become blocked in order to keep the GPU constantly working. If there aren’t enough threads available to the SMs, we just end up wasting performance.

Many of the other tests show similar patterns to the above as we increase the amount of work per-thread but often there is a sweet spot. Take a look at the results below, the fastest test is when we are processing 8 values per thread, anything above or below this is slower. Like with many things, it really is a delicate balancing act when you are trying to find the most optimal performance setup. We need to spawn enough threads to keep the GPU saturated, but the sweet spot here is where we can keep the GPU fully saturated but allow each thread to do as much work as possible.

Trace Metrics

512MiB vs 4MiB:

I hypothesised earlier that the performance difference I observed between the 512MiB and 4MiB tests were due to memory bandwidth, lets see if we can prove that using hardware metrics provided by Nvidia Nsight Systems.

We’ll compare the simple Add_AtomicPerThread<8,512> test with 512MiB x 1 run and 4MiB x 128 runs:

512MiB x 1

4MiB x 128

The results show as huge difference in the throughput between the two tests, with the 512MiB test at nearly 95% throughput on VRAM while the 4MiB test is closer to 1%. L1TEX and VidL2 throughput are both also much higher compared to the 512MiB, a good indication that the data we are processing is remaining in the caches and achieving better throughput as a result.

Add_SharedWithReduceAddSyncPerWarp vs Add_AtomicPerThread:

Next, lets compare two very different functions running the 4MiB x 128 test setup and see how the trace metrics look. On average, Add_SharedWithReduceAddSyncPerWarp is faster than Add_AtomicPerThread so it’s interesting to observe that Add_AtomicPerThread has better throughput for L1TEX as well as much better overall warp occupancy. I think it’s important to think about what the throughput stats actually mean, really it’s pointing to us the units which are under the most pressure. We saw above in the 512MiB vs 4MiB test that the VRAM throughput was 95% vs 1%, and this explains the huge difference in performance there, but just because we were achieving 95% throughput in that unit, does not mean it’s faster.

Add_SharedWithReduceAddSyncPerWarp< 4, 128 > (~1.7086ms)

Add_AtomicPerThread< 8, 512 > (~2.1138ms)

Conclusions

The main objective here was a personal objective of mine to introduce myself to CUDA programming, understand how it’s written, how to profile it and performance characteristics.

I think one thing that’s shown here is how in order to improve performance of your application, you must target what the primary bottleneck is. Simply reducing global atomic operations overall may not be the biggest improvement you can make, if it isn’t your primary bottleneck. Often, one of the most important things you can do is always try to keep the SM units fully saturated with as much work as possible and spread your workload evenly across as many threads as possible.

Something else that has stood out here is that using shared memory is very fast. Often the test scenarios which are more complex, involving more code and doing atomic writes to shared memory are amongst the fastest.

Further Investigations:

• Why is the occupancy so far off 100% for the traces above? What could be done to improve this further?
• Are there other CUDA features that could improve our parallel add further? Such as using CUDA Tiles or Clusters? Or Tensor features?
• How would these results differ on different GPUs?
• Are there more trace statistics we can gather to better understand our results, or better understand our existing trace data?

Source Code is accessible here: https://github.com/comeradealexi/CUDAParallelAddBenchmarks/blob/main/kernel.cu