A Worked CUDA Optimization: Debugging Tools and CUB Shift the Bottleneck from Median Kernels to Memory Allocation

NVIDIA CUDA remains the foundation of GPU-accelerated computing, powering everything from scientific simulations to large-scale AI training.

But writing correct, maintainable, and performant CUDA code can be challenging: memory bugs hide in plain sight, performance bottlenecks can be invisible without the right instrumentation, and hand-rolled GPU algorithms rarely match the efficiency of optimized libraries. Fortunately, the modern CUDA toolchain has matured significantly, and many of these challenges now have straightforward solutions.

In this blog post, we will walk through the tools NVIDIA offers to debug, benchmark, and improve your code. With only small line changes each time, we are going to make the example code safer, easier to maintain, and faster.

Across six incremental steps, this post will cover:

• How to easily find indexing bugs by adopting the modern CCCL API and Compute Sanitizer
• How to improve Nsight Systems benchmarks with NVTX
• How to use CUB’s optimized algorithms at the block and device level
• How to manage GPU memory through pooled containers
• How to speed up host-to-device transfers with pinned containers
• How to parallelize GPU work by giving each thread its own stream and asynchronous transfers

As a companion to this blog post, we provide the code and the option to run on Google Colab.

From an input stream of red, green, and blue images, start by transferring the data from the CPU to the GPU. Then convert those from RGB to grayscale.

Then, for each 32 by 32 pixel tile in the image, compute the median by sorting the pixels and selecting the middle value. Finally, copy the median of each tile back to the CPU.

Below is the full starting code. Each step in this post improves on it.

#define CUDA_CHECK_ERROR(call) do { \ cudaError_t err = call; \ if (err != cudaSuccess) { \ std::cerr << "CUDA error in " << __FILE__ << " at line " << __LINE__ << ": " \ << cudaGetErrorString(err) << std::endl; \ std::exit(EXIT_FAILURE); \ } \ } while (0) // Alias for an image pixel using pixel_t = uint8_t; // Kernel converting the red, green and blue images into a single gray image __global__ void computeRGBToGray(const pixel_t* d_image_r, const pixel_t* d_image_g, const pixel_t* d_image_b, pixel_t* d_image_gray, int width, int height) { // Compute the thread global index in the grid const int x = threadIdx.x + blockIdx.x * blockDim.x; const int y = threadIdx.y + blockIdx.y * blockDim.y; // Boundary check selecting only threads within the image boundary if (x < width && y < height) { // Compute the thread index in the image const int i = x + y * width; // Convert from rgb to grayscale and store the result in global memory d_image_gray[i] = static_cast<pixel_t>(0.299f * d_image_r[i] + 0.587f * d_image_g[i] + 0.114f * d_image_b[i]); } } // Kernel computing the median of each tile in the grayscale image template <int TILE_WIDTH, int HISTO_SIZE> __global__ void computeMedian(pixel_t *d_image_gray, pixel_t *d_median, int width, int height) { // Compute the thread global index in the grid const int x = threadIdx.x + blockIdx.x * blockDim.x; const int y = threadIdx.y + blockIdx.y * blockDim.y; // Boundary check selecting only threads within the image boundary if (!(x < width && y < height)) return; // Allocate the shared memory in which we will store the tile __shared__ pixel_t tile[TILE_WIDTH * TILE_WIDTH]; // Compute the thread index in the image const int index = x + y * width; // Load the tile's grayscale value from global memory into shared memory tile[index] = d_image_gray[index]; // Synchronize to make sure all threads have loaded their data __syncthreads(); // Sort the tile array using a single threaded bubble sort if (threadIdx.x == 0 && threadIdx.y == 0) { for (int i = 0; i < TILE_WIDTH * TILE_WIDTH; ++i) for (int j = i + 1; j < TILE_WIDTH * TILE_WIDTH; ++j) if (tile[i] > tile[j]) cuda::std::swap(tile[i], tile[j]); // Each thread block stores the median, found in the middle index after sorting, in the global median array const int medianIndex = (TILE_WIDTH * TILE_WIDTH) / 2; d_median[blockIdx.x + blockIdx.y * gridDim.x] = tile[medianIndex]; } } int main() { // Define all the example constants constexpr auto TILE_WIDTH = 32; constexpr auto HISTO_SIZE = 256; constexpr auto NB_TILE_X = 250; constexpr auto NB_TILE_Y = NB_TILE_X; constexpr auto IMAGE_LENGTH = TILE_WIDTH * NB_TILE_X; constexpr auto IMAGE_SIZE = IMAGE_LENGTH * IMAGE_LENGTH; constexpr auto NB_IMAGES = 3; constexpr auto INIT_VALUE = 4; // Allocate the CPU memory to store the images tiles medians and for the red, green, blue and grayscale images std::vector<std::vector<pixel_t>> h_images_r(NB_IMAGES, std::vector<pixel_t>(IMAGE_SIZE, 4)); std::vector<std::vector<pixel_t>> h_images_g(NB_IMAGES, std::vector<pixel_t>(IMAGE_SIZE, 4)); std::vector<std::vector<pixel_t>> h_images_b(NB_IMAGES, std::vector<pixel_t>(IMAGE_SIZE, 4)); std::vector<std::vector<pixel_t>> h_images_gray(NB_IMAGES, std::vector<pixel_t>(IMAGE_SIZE, 0)); std::vector<std::vector<pixel_t>> h_medians(NB_IMAGES, std::vector<pixel_t>(NB_TILE_X * NB_TILE_Y)); // Run the image processing pipeline for each image, in parallel #pragma omp parallel for for (int i = 0; i < NB_IMAGES; ++i) { pixel_t *d_image_r, *d_image_g, *d_image_b, *d_image_gray, *d_median; // Allocate the GPU memory for each container CUDA_CHECK_ERROR(cudaMalloc(&d_image_r, IMAGE_SIZE * sizeof(pixel_t))); CUDA_CHECK_ERROR(cudaMalloc(&d_image_g, IMAGE_SIZE * sizeof(pixel_t))); CUDA_CHECK_ERROR(cudaMalloc(&d_image_b, IMAGE_SIZE * sizeof(pixel_t))); CUDA_CHECK_ERROR(cudaMalloc(&d_image_gray, IMAGE_SIZE * sizeof(pixel_t))); CUDA_CHECK_ERROR(cudaMalloc(&d_median, (NB_TILE_X * NB_TILE_Y) * sizeof(pixel_t))); // Copy the memory of each container from CPU to GPU CUDA_CHECK_ERROR(cudaMemcpy(d_image_r, h_images_r[i].data(), IMAGE_SIZE * sizeof(pixel_t), cudaMemcpyHostToDevice)); CUDA_CHECK_ERROR(cudaMemcpy(d_image_g, h_images_g[i].data(), IMAGE_SIZE * sizeof(pixel_t), cudaMemcpyHostToDevice)); CUDA_CHECK_ERROR(cudaMemcpy(d_image_b, h_images_b[i].data(), IMAGE_SIZE * sizeof(pixel_t), cudaMemcpyHostToDevice)); // Launch a GPU kernel to convert the RGB images to grayscale dim3 blockSize(TILE_WIDTH, TILE_WIDTH); dim3 gridSize(cuda::ceil_div(IMAGE_LENGTH, blockSize.x), cuda::ceil_div(IMAGE_LENGTH, blockSize.y)); computeRGBToGray<<<gridSize, blockSize>>>(d_image_r, d_image_g, d_image_b, d_image_gray, IMAGE_LENGTH, IMAGE_LENGTH); CUDA_CHECK_ERROR(cudaGetLastError()); // Launch the GPU kernel to compute the median of every tile in the image computeMedian<TILE_WIDTH, HISTO_SIZE><<<gridSize, blockSize>>>(d_image_gray, d_median, IMAGE_LENGTH, IMAGE_LENGTH); CUDA_CHECK_ERROR(cudaGetLastError()); // Copy the GPU median memory back to the CPU CUDA_CHECK_ERROR(cudaMemcpy(h_medians[i].data(), d_median, (NB_TILE_X * NB_TILE_Y) * sizeof(pixel_t), cudaMemcpyDeviceToHost)); // Free the GPU memory CUDA_CHECK_ERROR(cudaFree(d_image_r)); CUDA_CHECK_ERROR(cudaFree(d_image_g)); CUDA_CHECK_ERROR(cudaFree(d_image_b)); CUDA_CHECK_ERROR(cudaFree(d_image_gray)); CUDA_CHECK_ERROR(cudaFree(d_median)); } return 0; }

This code starts by defining two kernels:

• computeRGBToGray loads the values of the red, green, and blue input images to convert and write them to the grayscale output image.
• computeMedian computes the median of each tile of an input grayscale image. Each thread block loads the tile from global memory to shared memory. Then, a single thread is used to sort the array and write the value found at the middle index, corresponding to the median, in the global output median array.

In the main, after defining the constants used for the example, the CPU memory is allocated for each image and for the medians.

The image processing pipeline is then run for each of the three images, in parallel, using OpenMP. The pipeline starts by allocating the required memory on the GPU before transferring the data from the CPU to the GPU. The two kernels to convert from RGB to grayscale and to compute the median are launched afterward. Finally, we copy back the median results to the CPU before freeing the memory.

This code has several flaws that will be addressed, step by step.

Let’s start by running the code.

code_steps$ ./build/0_base_error_example CUDA error in 0_base_error_example.cu at line 105: an illegal memory access was encountered

While the code has some error checking in it, when you get an error message like this “illegal memory access,” you should start by using compute-sanitizer to investigate further.

Using Compute Sanitizer, the NVIDIA functional correctness checking suite, we can directly identify a bug that can be hard to spot:

$ compute-sanitizer ./build/0_base_error_example ========= COMPUTE-SANITIZER ========= Invalid __shared__ write of size 1 bytes ========= at void computeMedian<(int)32, (int)256>(unsigned char *, unsigned char *, int, int)+0x170 in 0_base_error_example.cu:55 ========= by thread (0,3,0) in block (20,0,0) ========= Access at 0x6440 is out of bounds

Running the above directly shows that the code suffers from an out-of-bound shared write at line 55 of 0_base_error_example.cu .

tile[index] = d_image_gray[index];

The line incorrectly loads data in shared memory, using a global index. Since shared memory is defined at the thread block level, we need to change the indexing. To avoid indexing mistakes, a new API was introduced in CCCL to distinguish global and block-level indexing. To use it, you first need to launch your kernel using the new cuda::launch API:

auto config = cuda::make_config(cuda::block_dims(...), cuda::grid_dims(...)); cuda::launch(stream, config, kernel_name<decltype(config)>, input)

Then, using the new indexing API inside the kernel:

template <typename Configuration> __global__ void kernel_name(Configuration config, ...) { // Retrieve and expand each global index const auto [x, y, z] = cuda::gpu_thread.index(cuda::grid, config); // Retrieve the block index structure (containing block_idx.x, .y, .z) const auto block_idx = cuda::gpu_thread.index(cuda::block, config); }

Without using compute-sanitizer or the new API, this mistake could also have been directly spotted by using cuda::std::span or its n-dimensional variant cuda::std::mdspan instead of raw pointers. cuda::std::span and cuda::std::mdspan are non-owning views over contiguous memory and are useful to abstract the exact container away. Accessing data through spans is safer than through raw pointers in part because in debug mode, out-of-bounds accesses will trigger an assertion.

The kernel should be updated as:

// Alias for a 2-dimensional mdspan template <typename T> using span_2d = cuda::std::mdspan<T, cuda::std::dims<2>>; template <typename Configuration> __global__ void computeMedian(..., span_2d<const pixel_t> d_image_gray, ...)

If you run the code with these changes you’ll get something like the following:

$ ./build/1_span libcudacxx/include/cuda/std/__mdspan/mdspan.h:436: operator(): block: [16,0,0], thread: [0,30,0] Assertion `mdspan: operator() out of bounds access` failed.

Memory accesses in shared memory should also be protected by using a cuda::shared_memory_mdspan as in the following snippet:

__shared__ pixel_t shared[TILE_WIDTH * TILE_WIDTH]; cuda::shared_memory_mdspan tile_2d(shared, TILE_WIDTH, TILE_WIDTH);

Now you can run and everything should execute properly without errors.

Using the new launch API and its indexing mechanism, spans over raw pointers, and compute-sanitizer, out-of-bounds accesses either do not happen or are caught right away. For more information on compute-sanitizer, see Efficient CUDA Debugging: How to Hunt Bugs with NVIDIA Compute Sanitizer.

Now the code is bug-free, it is ready to be benchmarked using NVIDIA Nsight Systems. It allows you to visualize the program timeline: know when each function is called and for how long.

To make the timeline visualization easier, we wrap every interesting code section using NVTX:

void image_compute(...) { // NVTX range for the scope of the whole function nvtx3::scoped_range fun_scope("Image compute"); // NVTX range that is pushed and then popped for a specific code section nvtxRangePushA("Kernel median"); // Launch the GPU kernel to compute the median of every tile in the image ... // Pop the range at the end of the specific code section nvtxRangePop(); }

Which yields the following result:

In the GPU hardware (CUDA HW) section of the profiler output in Figure 3, above, it is reported that the GPU is mainly busy with kernels (98.5% of the GPU time) while the memory operations only take 1.5% of the GPU time.

Of the two kernels, the one computing the medians is taking the majority of the runtime with 2.1 seconds for each image (see the yellow box on the right, where the stats for computeMedian are shown, and the elapsed time is 2.142s).

From the CPU (thread) section, we see the image computation takes 6.8 seconds in total, with most of the time being spent on computing the medians for the three grayscale images.

We now know the first operation to optimize in order to have the greatest impact. For more information on Nsight Systems, see Optimizing CUDA Memory Transfers with NVIDIA Nsight Systems. For more information on NVTX, see CUDA Pro Tip: Generate Custom Application Profile Timelines with NVTX.

When dealing with common algorithms, writing custom kernels is error-prone and will most likely result in an inefficient implementation. Whenever possible, both for device-side patterns and for in-kernel primitives, it is recommended to use CUB.

CUB is the NVIDIA parallel algorithm library shipped through CCCL. It exposes highly optimized routines at several granularities: device-wide ( cub::Device* ), block-level ( cub::Block* ), and warp-level ( cub::Warp* ).

For the RGB to grayscale step, we can replace the custom kernel with cub::DeviceTransform::Transform . It applies a user-provided function to a tuple of input iterators and writes the result to an output iterator, executing on the GPU:

// Use CUB to convert the RGB images to grayscale cub::DeviceTransform::Transform( cuda::std::make_tuple(d_image_r, d_image_g, d_image_b), // inputs d_image_gray, // output IMAGE_SIZE, // size [] __host__ __device__ (pixel_t r, pixel_t g, pixel_t b) // functor { return static_cast<pixel_t>(0.299f * r + 0.587f * g + 0.114f * b); }, stream);

For the median, programming a parallel block-level sort by hand is complex and slow. Instead, we directly leverage CUB’s block-level radix sort inside the kernel:

// Declare and allocate the storage for CUB BlockRadixSort using BlockRadixSort = cub::BlockRadixSort<...>; __shared__ typename BlockRadixSort::TempStorage temp_storage; // Load the tile's grayscale value from global memory pixel_t thread_keys[1]; thread_keys[0] = d_image_gray(y, x); // Perform the thread-block-level radix sort BlockRadixSort(temp_storage).Sort(thread_keys); // Select the thread found at the middle index // Write its value which is, after sorting, the median, in the global median array if (block_idx.x == TILE_WIDTH / 2 && block_idx.y == TILE_WIDTH / 2) d_median(grid_block_idx.y, grid_block_idx.x) = thread_keys[0];

Following this change, we benchmark again using Nsight Systems:

The time required to compute the median is now only 773 microseconds (again, look at the elapsed time in the computeMedian yellow pop-out image), 2717x faster. The overall time to compute all three images is now 635 milliseconds, 10x faster.

If we reassess the current bottleneck: the time spent on memory allocations represents around 83% of the total image compute runtime.

This can be greatly improved.

Allocating GPU memory using cudaMalloc can have unexpected negative effects: leaks by forgetting to call cudaFree and costly memory operations in critical parts of your code.

Instead, we recommend using CCCL’s asynchronous memory containers, cuda::device_buffer . Like C++ std::vector , the memory is automatically deallocated once the container goes out of scope.

Additionally, a memory pool backs the buffer, so repeated allocations and deallocations do not pay the full cost of cudaMalloc / cudaFree every time.

To use the GPU memory containers, we update the code accordingly:

// Resource to handle the GPU memory allocations cuda::device_memory_pool_ref device_resource = cuda::device_default_memory_pool(cuda::device_ref{0}); // Explained at a later stage, unimportant for now cuda::stream stream{cuda::device_ref{0}}; // Allocate the GPU memory using uninitialized containers cuda::device_buffer<pixel_t> d_image_r = cuda::make_buffer<pixel_t>(stream, device_resource, IMAGE_SIZE, cuda::no_init); ...

After this change, we analyze the timeline again:

The time spent on memory allocation is now almost nonexistent; the time it takes to compute an image has improved by 2.6x.

The GPU time is now memory-dominated. Almost all the time to compute all images is spent on first copying the red, green and blue, for the three images, from CPU to GPU.

It is possible to greatly speed up the host-to-device memory transfers.

CPU data allocations are pageable by default, which the GPU cannot access directly. The CUDA driver must first allocate a temporary page-locked, or pinned, host array, copy the host data to the pinned array, and then transfer the data from the pinned array to device memory.

When it is known in advance that CPU memory will be copied to the GPU, it is advised to directly allocate using pinned memory.

CCCL exposes a pinned-memory host container, cuda::host_buffer , that can be built through the cuda::make_pinned_buffer factory:

// Allocate the CPU memory to store the image tiles, medians, and for the red, green, blue, and grayscale images // Those CPU containers, contrary to std::vector, are allocated using pinned memory std::vector<cuda::host_buffer<pixel_t>> h_images_r(NB_IMAGES, cuda::make_pinned_buffer<pixel_t>(stream, IMAGE_SIZE, ...)); ...

Following those changes, we can benchmark again:

The time it takes to do the host-to-device memory transfers has been significantly reduced; it now only takes 25 ms to compute all images, 10x faster.

For more information on pinned memory, see How to Optimize Data Transfers in CUDA C/C++.

One surprising behavior might have caught the eye of some readers since the beginning:

Though we are using different CPU threads, all operations (memory and kernels) are being executed sequentially on the GPU.

Let’s fix it.

By default, all operations (kernels, memory allocations, or transfers) are launched on what we call the default stream: it can be viewed as a queue of tasks the GPU needs to execute in order.

In this example, we need a stream for each image/thread. CCCL provides cuda::stream , an owning self-managed version of CUDA streams. It can simply be constructed inside the parallel for loop, so each OpenMP thread gets its own queue of GPU work.

To efficiently leverage streams, we also need to use the asynchronous API: each GPU operation launched by the CPU should not be waited upon until completion. To saturate the GPU, each CPU thread should launch as many operations as possible, as fast as possible, without waiting for them to first complete. Kernels and CUB device calls are already asynchronous by default and are launched on the passed stream. To launch asynchronous copies between the host and the device, we use the new CCCL cuda::copy_bytes API.

It is advised in any modern CUDA code never to rely on the default stream and to always rely on streams.

We update the code accordingly:

A dedicated init_stream is used for the initial allocations of the pinned host buffers. Each iteration of the parallel for loop now owns its own cuda::stream for the computation pipeline:

// Stream used for initial host buffer allocations cuda::stream init_stream{cuda::device_ref{0}}; ... // Allocate the host pinned buffers on init_stream: std::vector<cuda::host_buffer<pixel_t>> h_images_r(NB_IMAGES, cuda::make_pinned_buffer<pixel_t>(init_stream, IMAGE_SIZE, ...)); ... // Sync before launching operations on another stream: init_stream.sync(); #pragma omp parallel for for (int i = 0; i < NB_IMAGES; ++i) { // One different stream per thread cuda::stream stream{cuda::device_ref{0}}; ... // GPU buffer allocations using the stream owned by each thread cuda::device_buffer<pixel_t> d_image_r = cuda::make_buffer<pixel_t>(stream, device_resource, IMAGE_SIZE, cuda::no_init); ... // Copy the memory of each container from CPU to GPU asynchronously using the stream owned by each thread cuda::copy_bytes(stream, h_images_r[i], d_image_r); ... // Use CUB to convert the RGB images to grayscale asynchronously using the per thread stream cub::DeviceTransform::Transform(..., stream.get()); // Launch the GPU kernel to compute the median of every tile in the image using the per thread stream cuda::launch(stream, ...); // Copy the GPU median memory back to the CPU cuda::copy_bytes(stream, d_median, h_medians[i]); ... // To make sure the copy bytes is finished before accessing results on the host stream.sync(); }

Following those changes, we can take a final look at the timeline:

We now have a complete overlap between our kernels and memory copies.

The final duration to compute all three images following all our improvements is 23 milliseconds, starting from 6.8 seconds.

Using the CUDA Developer’s Toolbox, we made the code safer, easier to maintain and faster. No low-level optimizations were used, yet the code is 300x faster.

Try out this code yourself, and run on Google Colab if you like.

We have also built a full class to learn how to use those tools in detail. It is freely available on YouTube alongside links to practice on Google Colab.