CUDA Kernel Explained: Thread Hierarchy, Execution, and When to Write Your Own

A CUDA kernel is a C++ function marked __global__ that the CPU launches and the GPU runs many times in parallel, once per thread, where an ordinary C++ function runs once. You set the number of threads at launch with the triple-chevron syntax <<<blocks, threadsPerBlock>>>. Inside the kernel, each thread works out its own global index from blockIdx, blockDim and threadIdx, and that index tells it which data element it owns. Before the kernel can read anything, the data has to be in device memory, because the host and the device keep separate memory spaces. Writing your own kernel is warranted when no library covers the operation you need. For standard operations, NVIDIA’s tuned libraries are the default.

The rest of this article builds that answer from one working example. We start with the kernel itself, then cover how the hardware runs it, how to launch it, how data gets there, where it goes in a real project, and how to decide whether to write one.

A complete CUDA kernel example: vector addition

The example here is vector addition. NVIDIA’s own programming guide uses it to illustrate the built-in threadIdx variable. Here is the device side:

__global__ void vecAdd(const float* A, const float* B, float* C, int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // this thread's element
    if (i < n) {                                      // guard the tail
        C[i] = A[i] + B[i];
    }
}

There are three things to notice. The __global__ specifier marks the function as a kernel entry point: it runs on the device and can be called from the host. A __global__ function must return void and cannot be a member of a class. You write the code from the point of view of one thread, and the parallelism comes from launching enough threads.

Helper functions marked __device__ are compiled for the GPU and can be called from other __device__ or __global__ functions. NVIDIA’s Best Practices Guide suggests building a kernel from many short __device__ functions instead of one monolithic __global__ function, so each piece can be tested on its own.

How does each thread work out which element it handles?

Each thread combines three built-in variables into one global index: blockIdx.x * blockDim.x + threadIdx.x. threadIdx is the thread’s position within its block, blockIdx is the block’s position within the grid, and blockDim is the block size you chose at launch. So thread 5 of block 2, with 256 threads per block, handles element 2 × 256 + 5 = 517.

The if (i < n) check is not optional. You launch a whole number of blocks, so the grid usually holds more threads than there are elements. Without the guard, the threads past the end read and write memory that does not belong to the array.

The variables are three-component vectors, so the same scheme works in two or three dimensions. That fits images, matrices and volumes. A 2D kernel computes row = blockIdx.y * blockDim.y + threadIdx.y and col = blockIdx.x * blockDim.x + threadIdx.x, as in the programming guide’s matrix multiply, where each thread computes one element of C.

The mapping from index to address also affects memory performance. The hardware merges the global memory accesses of the threads in a warp into as few 32-byte transactions as it can. When consecutive threads touch consecutive elements, as they do in the vector add above, the accesses coalesce. NVIDIA calls proper coalescing one of the most important performance considerations for CUDA kernels. In a 2D block, threadIdx.x varies fastest, so it should be the index that walks along contiguous memory.

Thread hierarchy: how the GPU executes a launched kernel

A launch from the host creates one grid. The grid is made of thread blocks, and each block is made of threads. For execution, the hardware groups the threads of a block into warps of 32, which is the smallest unit of parallelism the device schedules. Once the kernel is launched, the scheduler assigns blocks to streaming multiprocessors (SMs). Your code cannot control or query which block runs where or in what order.

That freedom comes with a rule: blocks must be able to execute independently, in any order, in parallel or in series. The same rule is what lets one compiled program scale across GPUs with different numbers of multiprocessors. The practical consequence is that threads in the same block can synchronise and threads in different blocks cannot:

Within a block, threads can share data through __shared__ memory and wait for each other at a __syncthreads() barrier.
Across blocks, there is no grid-wide barrier in the basic model. Global-memory atomics give controlled read-modify-write access, and NVIDIA advises using them sparingly.
When threads in different blocks must share data, it goes through global memory in two kernel launches, one that writes and one that reads. That costs extra launches and extra memory traffic.

Branching inside a warp is a performance concern, not a correctness one. The programming guide says you can mostly ignore SIMT behaviour when you are designing for correctness, and that performance improves when threads in a warp seldom diverge. Advice like “avoid branches” only makes sense once a profile shows that divergence is what this kernel is actually paying for.

How do you launch a CUDA kernel with «<blocksPerGrid, threadsPerBlock»>?

You put the execution configuration in triple chevrons between the kernel name and its arguments. The first value is the grid size in blocks and the second is the block size in threads. The block count is a ceiling division of the problem size by the block size, which is what NVIDIA’s example computes with cuda::ceil_div:

int threads = 256;
int blocks  = (n + threads - 1) / threads;   // ceil(n / threads)
vecAdd<<<blocks, threads>>>(dA, dB, dC, n);

The full form is <<<Dg, Db, Ns, S>>>. Ns is the number of bytes of dynamically allocated shared memory per block and S is the stream. Both are optional. With 1D grids and blocks, plain integers work. For 2D or 3D layouts you pass dim3 values.

A block can hold at most 1024 threads on current GPUs, because all of a block’s threads are expected to reside on one SM and share its resources. The launch fails if the grid or block dimensions exceed the device maximums, or if the requested shared memory exceeds what is available.

How many threads per block, and why the maximum is often not the fastest

NVIDIA’s Best Practices Guide gives two starting points. Between 128 and 256 is a good initial range to experiment with. The grid should also contain more blocks than the GPU has multiprocessors, so that every SM gets work.

The maximum is not automatically best, because how many blocks fit on an SM at once depends on registers and shared memory as well as threads. NVIDIA’s programming guide gives an example for a device with 2048 resident threads per SM: 768-thread blocks allow only two blocks per SM, which caps occupancy at 75%. Register use also limits residency. In the guide’s compute capability 6.x example, one extra register per thread drops an SM from two 512-thread blocks to one. Occupancy, and how to tell whether it is your actual limiter, belongs to profiling GPU kernels; here it is enough to know that block size is a tuning parameter, not a constant.

Why a kernel fails silently, and how to make it speak

A kernel launch is asynchronous: control returns to the CPU thread before the kernel finishes, and possibly before it starts. A __global__ function must return void, and the call returns before the device has finished. Errors therefore show up in two places. A bad configuration is reported by cudaGetLastError() straight after the launch. A fault during execution is reported by the next synchronising call, such as cudaDeviceSynchronize(), which blocks until all earlier GPU work has completed. Every CUDA API call returns a cudaError_t, and NVIDIA’s guidance for production code is to check every one.

#define CUDA_CHECK(call)                                              \
    do {                                                              \
        cudaError_t e = (call);                                       \
        if (e != cudaSuccess) {                                       \
            fprintf(stderr, "%s:%d %s\n", __FILE__, __LINE__,         \
                    cudaGetErrorString(e));                           \
            exit(1);                                                  \
        }                                                             \
    } while (0)

vecAdd<<<blocks, threads>>>(dA, dB, dC, n);
CUDA_CHECK(cudaGetLastError());        // launch/configuration errors
CUDA_CHECK(cudaDeviceSynchronize());   // errors raised while running

Asynchrony also catches out naive timing. A CPU timer stopped right after the launch line measures the launch, not the kernel.

How do you allocate GPU memory and copy data around a kernel call?

The explicit sequence has five steps: allocate device memory with cudaMalloc, copy inputs over with cudaMemcpy(..., cudaMemcpyHostToDevice), launch the kernel, copy results back with cudaMemcpyDeviceToHost, and release the memory with cudaFree. Here it is around the vector add:

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

int main()
{
    const int n = 1 << 20;
    const size_t bytes = n * sizeof(float);
    std::vector<float> hA(n, 1.0f), hB(n, 2.0f), hC(n);

    float *dA, *dB, *dC;
    CUDA_CHECK(cudaMalloc(&dA, bytes));
    CUDA_CHECK(cudaMalloc(&dB, bytes));
    CUDA_CHECK(cudaMalloc(&dC, bytes));

    CUDA_CHECK(cudaMemcpy(dA, hA.data(), bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(dB, hB.data(), bytes, cudaMemcpyHostToDevice));

    int threads = 256;
    int blocks  = (n + threads - 1) / threads;
    vecAdd<<<blocks, threads>>>(dA, dB, dC, n);
    CUDA_CHECK(cudaGetLastError());

    CUDA_CHECK(cudaMemcpy(hC.data(), dC, bytes, cudaMemcpyDeviceToHost));

    CUDA_CHECK(cudaFree(dA));
    CUDA_CHECK(cudaFree(dB));
    CUDA_CHECK(cudaFree(dC));
    printf("C[0] = %f\n", hC[0]);
}

The final copy does not need its own synchronisation call. cudaMemcpy is blocking, and work issued to the default stream runs in order, so the copy starts only after the kernel has finished.

The alternative is unified memory. cudaMallocManaged returns one pointer that both CPU and GPU code can use, and the driver moves the data on demand. It produces shorter code than explicit management. NVIDIA presents explicit management as the trade the other way: more verbose, but with control over when copies happen and where memory lives, and with the option of overlapping transfers with computation. With unified memory, cudaMemPrefetchAsync can move data to the GPU before the kernel needs it. If you stay explicit, page-locked host buffers from cudaMallocHost improve copy performance and are required for asynchronous transfers. Page-lock only the buffers you actually transfer, because pinning too much host memory can degrade performance on some systems.

Whichever style you use, host-device transfers are costly and should be minimised. Host-device copies have much lower bandwidth than device memory itself. NVIDIA’s guidance is to keep data on the device between kernel calls, even if that means running a step on the GPU that the CPU could do faster. Whether that cost pays off against a CPU implementation is the subject of when GPU offload actually beats a CPU implementation.

Where should CUDA kernels be defined in a C++ program with multiple source files?

Put kernels in .cu files. nvcc treats .cu as the extension for files that contain device code, and it accepts .cpp, .cc, .cxx and .c for host-only code. nvcc separates device code from host code and sends each to the right compiler. Host-only files can be built with your normal host compiler and linked with nvcc’s objects afterwards. Device-side declarations shared between .cu files conventionally go in .cuh headers.

A layout that keeps the rest of the project free of CUDA syntax works well. Define the kernel and a plain C++ launcher function in a .cu file, declare only the launcher in an ordinary header, and call it from .cpp code compiled by g++ or MSVC. The host compiler never sees a triple chevron.

The constraint appears when a kernel in one .cu file calls a __device__ function defined in another. By default nvcc uses whole-program compilation, which expects all GPU code a kernel uses to be in its own compilation unit. To link device code across files you compile with -dc (an alias for -rdc=true), which nvcc calls separate compilation. Your CPU compiler cannot link CUDA device code, so the build needs a step in which nvcc links it:

nvcc -dc  kernels.cu helpers.cu          # relocatable device code
nvcc -dlink kernels.o helpers.o -o gpu_link.o
g++ main.cpp kernels.o helpers.o gpu_link.o -lcudart -o app

Separate compilation can improve compile times and produce smaller binaries. It has a cost: when a call site and its callee sit in different compilation units, the compiler must follow the ABI and cannot optimise across the call. NVIDIA says the performance effect varies. For that reason we would keep small, hot __device__ helpers in a header included by the kernel’s own .cu file, so they can be inlined, and use cross-file device linking for code organisation rather than for inner-loop functions.

When is it worth writing your own CUDA kernel instead of using a library?

Default to a library and write a kernel when no library covers what you need. NVIDIA’s guidance points the same way across several documents. The programming guide tells developers to prefer the tuned CUDA Core Compute Libraries (CCCL) primitives over reimplementing algorithms such as reductions, for both productivity and performance. The Best Practices Guide names calling cuBLAS, cuFFT or Thrust as one way to parallelise a hotspot. For an advanced hardware feature set the guide describes, NVIDIA strongly recommends going through cuBLAS, cuDNN or cuFFT, or through CUTLASS inside device code.

Libraries win partly because a library is many kernels, not one. Simon Boehm’s SGEMM optimisation worklog notes that cuBLAS contains hundreds of SGEMM implementations and picks one at runtime from the matrix dimensions. He found 16 different kernels in use for square matrices up to 4096. After shared-memory tiling, his kernel reached about 2,200 GFLOPS, still far from the roughly 30 TFLOPS the GPU could deliver. The next kernel, which computes more results per thread, reached about 8,600 GFLOPS.

A custom kernel earns its place in a few specific situations:

No library implements the operation. NVIDIA names a missing library as the case where you may have to write your own GPU code in a high-level language such as C++.
The library cannot fuse what you need. NVIDIA’s handwritten-PTX article benchmarks a CUTLASS kernel that fuses GEMM with top-k and softmax in one pass.
A profile shows a hotspot that custom code would address. The Best Practices Guide frames a custom kernel as the step after you have located the hotspot and decided that custom code is the best approach.

Below CUDA C++ sits handwritten PTX, NVIDIA’s GPU assembly language, which NVIDIA describes as a last-resort optimisation. In its benchmark, handwritten PTX improved performance by 7% to 14% over the CUDA C++ version. NVIDIA stresses that the point is that gains are possible in carefully chosen situations, not the specific numbers. It also warns that low-level gains may not carry over to other GPU architectures, and that this kind of code adds development and debugging complexity.

Decision table: library, custom kernel or PTX

Situation	Reasonable default	Why
Standard dense maths (GEMM, FFT, convolution)	cuBLAS, cuFFT, cuDNN	Tuned, with many shape-specific kernels selected at runtime
Reductions and similar primitives	CCCL	NVIDIA recommends these over reimplementation
GEMM-like work with a custom epilogue or fused step	CUTLASS	Templates expose more control than cuBLAS, at the price of more code
Operation no library provides	Custom CUDA C++ kernel	NVIDIA’s stated case for writing your own
Chain of small ops round-tripping through global memory	Custom fused kernel, after profiling	Keeping data on-chip avoids memory traffic; worth it only if the profile shows the traffic dominates
Narrow, carefully chosen hotspots	Inline PTX with a C++ fallback	Measured gains in narrow cases, with portability risk

The common thread is measurement. NVIDIA’s programming guide says optimisation should be directed by measuring the performance limiters, because optimising instruction use in a memory-bound kernel gains nothing significant. Lower-level access opens more optimisation, but it pays off only with real commitment, and whether it is worth it is a profiling result, not a rule of thumb. Before you write custom code, the profiling-first approach to improving GPU performance is the better starting point than a hand-written kernel, and the Nsight Systems and Nsight Compute workflow is how that profile gets taken.

Your first kernel is a correctness baseline

The vector add above is correct, coalesced and checked for errors, and that is all it is. It does not show which performance limiter a real workload hits. Which optimisation pays off depends on the kernel’s performance limiters, and NVIDIA says to find them by measuring, for example with the CUDA profiler.

Before you rewrite anything, settle one question: is the time going to the kernel at all, or to the copies and launches around it? For short kernels, launch overhead can be a significant share of end-to-end time. Host-device copies can outweigh a cheap computation entirely. Answering that comes before any kernel tuning, and it is what profiling GPU kernels to find the real bottleneck walks through.