Part 2 — Basic Kernels

Everything in Part 1 ran on the host and called into a library. This part is about the code that runs on the GPU itself: how to write it, how to launch it, and how the threads that run it are arranged.

2.1 A first kernel

Let's start with a simple example: saxpy — single-precision a times x plus y, a basic BLAS operation. Here, the inputs are two floating-point arrays x and y of length n, as well as a single float a. The values of y are updated by adding scaled versions of x to them, such that y[i] <- a * x[i] + y[i].

This problem is specifically chosen to be trivially parallelizable: we can write down the equation for a single data element, and then only need to apply the function — or kernel in CUDA terminology — to every position in the array.

Device code has to be compiled by nvcc, which conventionally means putting it in a file with the extension .cu instead of .cpp.¹ To tell the NVCC compiler that a function is to be compiled for execution on the GPU instead of the CPU, the function needs to be marked with global or device. A __global__ function defines a GPU kernel, an entry-point that can only be called by a kernel launch. On the other hand, a __device__ function cannot be directly called from host code, but can be called from other GPU code (__device__ or __global__ functions).

Thus, the equation above can be coded as

__device__ void saxpy_impl(int i, float a, const float* x, float* y) {
    y[i] = a * x[i] + y[i];
}

Work partitioning and kernel launch

Now all that is missing is that we assign the different indices i to different worker threads on the GPU. In contrast to a CPU, where a program might run tens to at most hundreds of threads, and where thread creation is a heavyweight operation, a GPU can easily run tens of thousands of threads in parallel. In the next part, we will look more closely at how the GPU can manage that efficiently; for now, just assume that running this many threads is not only possible, but intended.

In this case, the simplest possible work-partitioning that we can do is to assign one thread to every index i. Mapping the identity of a thread to a unit of work is such a fundamental and frequent pattern of GPU programming that CUDA provides dedicated built-in variables for this.

The variable that carries that identity is threadIdx.x (.y and .z exist for multidimensional kernels), and reading it is all the first version of the kernel needs:

__global__ void saxpy_kernel(float a, const float* x, float* y) {
    int i = threadIdx.x;
    saxpy_impl(i, a, x, y);
}

Now, the only thing missing is how to launch the kernel and tell it how many threads to use. While there exists a function in the CUDA runtime API that can do this, most of the time we will use the convenient shorthand syntax <<<>>>, which NVCC automatically translates to the right API functions.

__host__ void saxpy_launcher(int n, float a, const float* x, float* y) {
    saxpy_kernel<<<1, n>>>(a, x, y);
}

This is now a fully-functioning kernel, but it only works for n <= 1024, using a tiny fraction of the GPU's compute power. To extend it to larger sizes, we need to consider one more detail about how the GPU executes a kernel.

2.2 Threads, Blocks, and Grids

Back in §1.4 The memory hierarchy, we showed that the GPU is partitioned into many streaming multiprocessors (SM). Threads follow that organization, being grouped together into Cooperative Thread Arrays (CTAs), or blocks.²

There are good reasons for this hierarchical organization of threads, as opposed to having just a flat index. It allows us to make guarantees about the subset of threads within one block that cannot be made for the full grid, i.e., the set of all threads created by the kernel launch. Threads within one block are tightly coupled both temporally and spatially:

• All threads within one block are executed at the same time.
• All threads within one block are executed on the same SM.

Together, these two properties imply that threads of a block can synchronize and communicate with each other, including the use of shared memory as a scratchpad. On the flipside, these guarantees put a limit on how many threads a block can hold until the SM runs out of resources. This limit depends on the amount of resources consumed by each thread, but on most GPUs there is an upper limit of 1024 threads per block even for the most frugal of threads.

In contrast, there are no such guarantees for threads belonging to different blocks. It is entirely possible that one block of threads runs and finishes its work before another one is even started. In fact, because there is a limit of how many blocks can be active on the GPU at the same time, it might even be required that some blocks finish before others can start. Thus, any attempt to synchronize between such blocks would just halt the GPU forever in a deadlock.

Which block ends up on which SM, and in what order, is not specified and left to the discretion of the runtime scheduler.

Every thread runs the same kernel from the top, so the only thing distinguishing one from another is where it sits in this arrangement: which thread it is within its block, and which block that is within the grid. Both, together with the sizes of the two groupings, come from four built-in variables:

All four are three-dimensional, with members .x, .y and .z; the other two dimensions exist to make grids over two- and three-dimensional data easier to write, and are of no concern for this part.

With these variables, each thread can work out which data element to access as illustrated below:

2.3 Automatic scaling

A convenient consequence of having the grid split into blocks, with blocks not required to be all in existence at the same time, allows us to adapt the number of threads to the problem size, independently of the number of SMs the GPU has available.

In the example above, we launched 1024 threads for 1024 elements, using only a single SM while the rest of the GPU idles. If instead our input had 102400 elements, we could have launched the kernel as saxpy_kernel<<<100, 1024>>>(a, x, y); and produced 100 blocks. Assuming, for the sake of argument, that the GPU had 100 SMs, that would be exactly one block per SM. Doubling the input size to 204800 elements, the code still works. The GPU schedules the first 100 blocks to its SMs, waits for them to finish, and then launches the next 100 blocks. To process the entire kernel, we now had to launch two blocks sequentially; the kernel execution is said to be in two waves.

What if we switched to a smaller GPU with only 50 SMs? The same kernel launch still gives correct results, only now it needs four waves.

The cost of this flexibility in scale is that we lose the ability to specify exact grid sizes: The number of threads in the grid must be a multiple of the block size. If the problem size is not an exact match, that means that we need to handle partial blocks: Spawn sufficiently many blocks that the entire problem is covered, and inside the kernel implement a guard that ensures that threads outside of the valid range of elements do not access any data.

2.4 Launching a kernel

The general kernel launch thus looks like this:

saxpy_kernel<<<num_blocks, block_size>>>(a, x, y);

The first parameter is the size of the grid, the second the size of a block; these are exactly the gridDim and blockDim the threads will read, and the number of threads started is their product.

Choosing the block size. Anything from 1 to 1024, but for reasons we will discuss in the next part, it is recommended to pick a multiple of 32. Often, 128 and 256 are reasonable defaults.

Choosing the grid size. Enough blocks to cover the array, which means dividing by the block size and rounding up — in integer arithmetic:

int num_blocks = (n + block_size - 1) / block_size;

Arguments are copied to the device. The values in the parameter list are transferred at call time, so they have to be values the device can use: scalars are fine, and pointers must point into device memory. Passing a host pointer compiles without complaint and faults when a thread dereferences it. The parameters are copied by sending their bit-representation to the GPU, so anything with non-trivial copy semantics (e.g., a std::vector) is not supported.

The launch is asynchronous. <<<>>> returns as soon as the work has been handed to the driver and the parameters have been copied and the host carries on. This means you cannot measure a kernel's execution time by reading the host clock on either side of the launch. It also means that you cannot get feedback from the kernel's execution directly after the launch.

Checking for errors. A launch has no return value to check, so the runtime records failures for later collection, and there are two separate moments at which something can go wrong:

saxpy_kernel<<<num_blocks, block_size>>>(a, x, y);
CUDA_CHECK(cudaGetLastError());       // did the launch get accepted?
CUDA_CHECK(cudaDeviceSynchronize());  // wait, then: did the threads run cleanly?

cudaGetLastError returns immediately and reports what the runtime could tell at submission time — a block larger than 1024 threads, a kernel that was never compiled for this architecture. Faults committed by the threads themselves, such as an out-of-bounds access, may not exist yet at that point. You need to wait for the kernel execution to complete, as implemented by cudaDeviceSynchronize, before those are available.

export CUDA_LAUNCH_BLOCKING=1

2.5 What nvcc actually does

A .cu file holds two languages at once, and nvcc is mostly a driver program that separates them. Host code is handed to the ordinary system compiler, with launches replaced by CUDA runtime calls.

Device code is first compiled to PTX, a virtual instruction set: it is assembly, but for an idealised GPU with unlimited registers and no particular architecture. PTX is then assembled by ptxas into SASS, the real machine code for one specific GPU architecture. Later in the course, we will look at examples of both PTX and SASS.

ptxas writes the SASS into a cubin, and nvcc collects every cubin it was asked to build, together with the PTX they were assembled from, into one fatbin that is embedded in the executable next to the ordinary host code.

The split into PTX and SASS allows for much simplified code portability: On x86, modern CPUs still carry a large amount of legacy instructions to ensure old programs remain runnable. For the GPU, instead, the executable includes the PTX version of the code, and the CUDA driver can generate SASS on the fly, for the hardware it actually runs on.

If SASS is generated by the driver, then why does NVCC still include it in the binary? Two reasons: First, running ptxas for every kernel on its first invocation adds runtime overhead that can be avoided by having the SASS readily available. Second, while a large amount of optimizations are performed by ptxas, certain choices need to be made already at the PTX level, such as the decision to use certain hardware features which are only present in some GPU models. Thus, common practice is to include SASS for the GPUs you currently care about, and PTX to ensure the code remains portable.

Which GPUs those are is set with -arch. -arch=sm_90 compiles for a Hopper card: it puts an sm_90 cubin in the fatbin and the PTX beside it, so the binary runs at full speed on Hopper and still runs on hardware released after it. Every architecture has such a number — sm_80 for Ampere, sm_100 for datacentre Blackwell and sm_120 for the consumer RTX50xx cards — and if the machine you compile on is the machine you run on, -arch=native looks the number up from the card that is installed.

nvcc -gencode arch=compute_80,code=sm_80 \
     -gencode arch=compute_90,code=sm_90 \
     -gencode arch=compute_100,code=sm_100 kernel.cu

2.6 Debugging kernel errors

As noted above, the kernel we show only works correctly when the number of elements is a multiple of the block size. Fixing this is your task in the exercise below. Here, we demonstrate instead which tools to use to debug kernels. These become almost indispensable for larger kernels with less-obvious deficiencies.

For reference, this is the full program we consider in this section:

#include <cuda_runtime.h>
#include <cstdio>
#include <vector>

inline void cuda_check_impl(cudaError_t status, const char* statement,
                            const char* file, int line) {
    if (status != cudaSuccess) {
        fprintf(stderr, "CUDA error in %s:%d (%s): %s: %s\n", file, line,
                statement, cudaGetErrorName(status), cudaGetErrorString(status));
        exit(EXIT_FAILURE);
    }
}
#define CUDA_CHECK(statement) cuda_check_impl((statement), #statement, __FILE__, __LINE__)

__device__ void saxpy_impl(int i, float a, const float* x, float* y) {
    y[i] = a * x[i] + y[i];
}

__global__ void saxpy_kernel(float a, const float* x, float* y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    saxpy_impl(i, a, x, y);
}

void saxpy_launcher(int n, float a, const float* x, float* y) {
    constexpr int BLOCK = 256;
    int num_blocks = (n + BLOCK - 1) / BLOCK;
    saxpy_kernel<<<num_blocks, BLOCK>>>(a, x, y);
}

int main() {
    // 1023 is one short of four whole blocks of 256, so the last block runs
    // with a single thread past the end of the arrays.
    const int N     = 1023;
    const float A   = 2.0f;
    const size_t sz = N * sizeof(float);

    std::vector<float> hX(N), hY(N);
    for (int i = 0; i < N; i++) {
        hX[i] = 0.5f * i;
        hY[i] = 1.0f;
    }

    float *dX, *dY;
    CUDA_CHECK(cudaMalloc(&dX, sz));
    CUDA_CHECK(cudaMalloc(&dY, sz));
    CUDA_CHECK(cudaMemcpy(dX, hX.data(), sz, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(dY, hY.data(), sz, cudaMemcpyHostToDevice));

    saxpy_launcher(N, A, dX, dY);
    CUDA_CHECK(cudaGetLastError());      // errors from the launch itself
    CUDA_CHECK(cudaDeviceSynchronize()); // errors from kernel execution

    CUDA_CHECK(cudaMemcpy(hY.data(), dY, sz, cudaMemcpyDeviceToHost));
    for (int i : {0, 1, N - 1})
        printf("y[%4d] = %g\n", i, (double)hY[i]);

    CUDA_CHECK(cudaFree(dX));
    CUDA_CHECK(cudaFree(dY));
    return 0;
}

It runs, prints the right numbers, and exits zero. Both error checks pass. The thread past the end reads and writes four bytes beyond a 4092-byte allocation, so everything is fine? No, we are just getting lucky here. Memory access protections are only enforced at the granularity of full memory pages, which are of size 2 MiB on the GPU. Thus, our beyond-the-end access still lands within the allocated page, reading and writing garbage. Only if you're unlucky, and someone passes your kernel an input array that ends just at a page boundary, would this code crash.

Without any additional tools, such infrequent bugs are hard to spot and even harder to reproduce. For this reason, CUDA ships four mechanisms to aid in debugging. In increasing order of effort, these are:

• Environment variables. Setting CUDA_LAUNCH_BLOCKING=1 disables the asynchrony of kernel launches, making it much easier to localize which kernel was responsible for a crash. The driver log, switched on with the CUDA_LOG_FILE environment variable, records additional information for failing API calls.
• print. The cuda runtime provides a version of printf that can be called from device code.
• compute-sanitizer re-runs the program with the device's memory accesses checked, helping to detect and pinpoint memory errors.
• cuda-gdb puts breakpoints into device code and steps through it.

The driver log needs no rebuild and no change to the program. Setting CUDA_LOG_FILE to stdout, stderr, or a path makes the driver write one line per failure as it happens:

[13:29:05.238][139569331646464][CUDA][E] Returning 2 (CUDA_ERROR_OUT_OF_MEMORY) from cuMemAlloc_v2

It needs driver r570 or newer, and it reports what the driver saw, so it complements CUDA_CHECK rather than replacing it. On the program above it stays silent: no API call fails, and the bug is inside a kernel that the driver considers to have run to completion.

printing With its simplicity, printf is a time-honoured tool for inspecting the state of programs during debugging and development. There are two major pitfalls to be aware of with printing from kernel code, however:

1. It does not actually print directly, but instead writes to a device-side buffer that is only later flushed to the host, so a kernel crash might hide the result of the last print statements. Not seeing a print does not necessarily imply that the print statement was not reached!
2. An unguarded print statement will be executed by thousands of threads, spamming the standard output. Typically, you'll want to focus on only a small subset of threads, and guard the print call with a conditional.

In this example, we could print the offsets the kernel tries to access, and see that one of them is past the range of valid values. There is a much better approach for these kinds of memory errors, though.

compute-sanitizer. This tool runs the program with every device memory access checked, against the bounds of the allocation rather than those of the memory page. It detects the problematic access:

$ compute-sanitizer --tool memcheck ./saxpy
========= COMPUTE-SANITIZER
========= Invalid __global__ read of size 4 bytes
=========     at saxpy_impl(int, float, const float *, float *)+0xa0 in saxpy.cu:28
=========     by thread (255,0,0) in block (3,0,0)
=========     Access to 0x7ef90c800ffc is out of bounds
=========     and is 1 bytes after the nearest allocation at 0x7ef90c800000 of size 4.092 bytes
=========         Device Frame: saxpy_kernel(float, const float *, float *)+0x50 in saxpy.cu:33
...
========= ERROR SUMMARY: 2 errors

The access is a 4-byte read; the culprit is thread 255 of block 3, the one thread the arithmetic put past the end; the address is one byte beyond the allocation; and the code is at saxpy.cu:28, inlined into the kernel at line 33.

Checking costs execution time. A kernel under the sanitizer runs an order of magnitude slower. Memory checking is the default. Other tools are --tool racecheck, initcheck and synccheck, which look for shared-memory races, uninitialised reads and mismatched barriers.

The file and line come from -lineinfo. Without it, the same run reports only

=========     at saxpy_kernel(float, const float *, float *)+0xa0

No source location, and the inlined saxpy_impl frame folded into its caller. -lineinfo adds a line table, but contrary to full debug information -G, it does not change the generated code at all, the only cost is a larger binary. As such, it is a good default during development, especially as lineinfo is also useful for profiling.

cuda-gdb is a full-fledged debugger, an extension to gdb to include GPU threads. It lets you stop inside a kernel, inspect a variable and step. In order for it to have full debugging capabilities, it wants -G, the heavier relative of -lineinfo: full debug information, generated by turning device optimization off. Unfortunately, that means this will be a substantially different execution than the actual release kernel, so bugs that only manifest due to specific compiler optimizations can remain untriggered during debug.

Also, given the massive number of threads running concurrently, the stepping experience can be quite different than on a CPU program. If you're debugging an actual kernel crash, however, running under cuda-gdb will trap at the location of the crash, giving you a chance to inspect the state of the program while already pointing you to the offending thread.

void axpy(int n, float alpha, const float* x, float* y);