Understanding the Problem: Variable-Sized Sub-Arrays in CUDA

In CUDA programming, a common challenge involves handling ragged arrays or flat arrays containing multiple sub-arrays of non-equal sizes. Suppose you have a large flattened key array and value array along with an offset/starts array marking the boundary of each sub-array. Sorting each segment independently and in parallel across the GPU is a crucial step for many parallel algorithms like graph processing, sparse matrix operations, and clustering.

A naive attempt often involves writing a custom CUDA kernel where each thread calls thrust::sort_by_key:

__global__ void sort_arrays(int* arrays, float* keys, int* starts) {
    const unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    // ...
    thrust::sort_by_key(thrust::device, keys + start, keys + start + size, arrays + start);
}

Why Thrust Inside CUDA Kernels Doesn't Work as Expected

Calling thrust::sort_by_key inside a custom device kernel using thrust::device either triggers CUDA Dynamic Parallelism (CDP) or falls back to sequential device execution (thrust::seq). This leads to several major issues:

  • Kernel Launch Overhead: Nesting kernel launches inside device code via CUDA Dynamic Parallelism introduces massive execution overhead, especially for thousands of tiny sub-arrays.
  • Sequential Fallback: If Thrust executes sequentially per thread, you lose parallel speedup across individual elements within a segment.
  • Suboptimal Thread Utilization: Assigning one thread per sub-array results in warp divergence and underutilizes GPU hardware for large segments.

The Solution: CUB Segmented Sort

The standard, highly optimized way to perform parallel sorting on ragged arrays in CUDA is by using CUB (CUDA Unbound), which is included directly in the NVIDIA CUDA Toolkit and integrated with modern Thrust releases.

Specifically, cub::DeviceSegmentedSort::SortPairsByKey allows you to sort multiple contiguous segments of keys and values simultaneously in a single, batched operation without writing custom device-side loops.

Implementation Example Using CUB

Here is how to set up and execute segmented sorting using CUB:

#include <cub/cub.cuh>
#include <cuda_runtime.h>

void sort_subarrays(
    float* d_keys_in,
    float* d_keys_out,
    int* d_values_in,
    int* d_values_out,
    int num_items,
    int num_segments,
    int* d_offsets // Array of size num_segments + 1 containing segment start indices
) {
    // Determine temporary device storage requirements
    void* d_temp_storage = nullptr;
    size_t temp_storage_bytes = 0;

    // Query temporary storage requirement
    cub::DeviceSegmentedSort::SortPairsByKey(
        d_temp_storage, temp_storage_bytes,
        d_keys_in, d_keys_out,
        d_values_in, d_values_out,
        num_items, num_segments,
        d_offsets, d_offsets + 1
    );

    // Allocate temporary storage
    cudaMalloc(&d_temp_storage, temp_storage_bytes);

    // Execute segmented sort
    cub::DeviceSegmentedSort::SortPairsByKey(
        d_temp_storage, temp_storage_bytes,
        d_keys_in, d_keys_out,
        d_values_in, d_values_out,
        num_items, num_segments,
        d_offsets, d_offsets + 1
    );

    // Free temporary memory
    cudaFree(d_temp_storage);
}

Alternative Approaches

Depending on your segment sizes, here are two alternative methods:

  • CUDA Streams: If you have a small number of very large sub-arrays, you can assign each sub-array sort to a different cudaStream_t on the host side using host-side Thrust calls.
  • CUB Block/Warp Primitives: If all your segments are guaranteed to be extremely small (e.g., fitting within a warp or block of 32–1024 elements), block-level or warp-level CUB sort primitives inside a single custom kernel will deliver maximum performance.

Conclusion

Do not invoke high-level Thrust routines directly inside a GPU kernel when trying to sort multiple sub-arrays. Instead, leverage cub::DeviceSegmentedSort to achieve maximum GPU throughput for non-equally sized sub-arrays.