GPU kernel basics (CUDA)

A GPU runs one small function — a kernel — on tens of thousands of threads at once. Writing a correct kernel is mostly ordinary C++; writing a fast one, and one that stays correct, comes down to a handful of rules about how those threads are grouped, how they touch memory, and where they must wait for each other. Those rules are what rendering and ML interviews ask about ("why is this transpose slow?", "what's wrong with this reduction?"), and what you need in the first week of any job with a CUDA backend.

There is no GPU in this sandbox, and no nvcc. Real CUDA appears here as fragments (they are what you'd write on a machine with a GPU), and the mechanics are demonstrated with ordinary C++ that computes exactly what the hardware rules say: how many memory transactions a warp issues, which shared memory bank a lane hits, how a block reduction proceeds. The build task gives you a deterministic CPU kernel simulator that measures all of this; you'll run real GPU code in the project when hardware is available.

1. The execution model: grid, block, warp, thread

A launch kernel<<<grid, block>>>(args) creates a grid of blocks; every block has the same number of threads (up to 1024). Both are dim3s, so you can lay them out in 1, 2 or 3 dimensions — 2D for images is the common case.

• A block runs on one SM (streaming multiprocessor) from start to finish. Its threads can share fast on-chip shared memory and wait for each other with __syncthreads(). Different blocks can't synchronise with each other inside a kernel and may run in any order, concurrently or not.
• The hardware schedules a block in warps of 32 threads with consecutive linear thread index (threadIdx.x varies fastest, then .y, then .z). A warp issues one instruction for its 32 lanes at a time (SIMT). If lanes disagree on a branch, the warp runs both paths with the inactive lanes masked off — divergence costs time, not correctness.
• Inside the kernel, the built-in variables gridDim, blockDim, blockIdx and threadIdx tell a thread where it is.

2. Index math, bounds checks and grid-stride loops

The first line of almost every kernel computes the element the thread owns:

__global__ void scale(float* v, int n, float k) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // global index
    if (i < n) v[i] *= k;                             // the bounds check
}
// host: enough blocks to cover n (round up!)
int threads = 256;
int blocks = (n + threads - 1) / threads;
scale<<<blocks, threads>>>(d_v, n, 2.0f);

The grid is rounded up, so the last block usually has threads past the end: the bounds check is not optional. Rounding down silently skips the tail. In 2D, each coordinate gets the same treatment — and both must be checked before computing the linear index, because y * width + x with x >= width is a perfectly valid index into the next row:

__global__ void shade(float4* img, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;   // threadIdx.x → x: see §4
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= width || y >= height) return;
    img[y * width + x] = make_float4(x, y, 0, 1);
}
dim3 block(32, 8);
dim3 grid((width + 31) / 32, (height + 7) / 8);
shade<<<grid, block>>>(d_img, width, height);

A grid-stride loop decouples the grid from the data size: launch a fixed number of blocks (say, a few per SM) and let each thread step by the whole grid. Any n works, consecutive lanes still touch consecutive elements, and the same kernel can be run with one block while debugging.

__global__ void saxpy(int n, float a, const float* x, float* y) {
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += gridDim.x * blockDim.x)
        y[i] = a * x[i] + y[i];
}

Here is the same arithmetic on the CPU — which thread handles which element for a deliberately tiny launch:

#include <print>

int main() {
    const int n = 20, block_dim = 4, grid_dim = 2;   // 8 threads for 20 elements
    for (int b = 0; b < grid_dim; ++b)
        for (int t = 0; t < block_dim; ++t) {
            std::print("block {} thread {}:", b, t);
            for (int i = b * block_dim + t; i < n; i += grid_dim * block_dim)
                std::print(" {}", i);
            std::println();
        }
    // A one-element-per-thread launch instead rounds the grid up:
    const int threads = 256, elems = 1000;
    int blocks = (elems + threads - 1) / threads;
    std::println("{} elements, {} threads/block: {} blocks, {} idle threads",
                 elems, threads, blocks, blocks * threads - elems);
}

Use int (or unsigned) for indices only when the product fits: y * width + x for a 50 000 × 50 000 image overflows 32 bits. Large buffers want size_t (or int64_t) arithmetic.

3. The memory hierarchy

The game is to read global memory as few times as possible, in as few transactions as possible, and to do reuse in registers and shared memory.

4. Coalescing: what a warp's 32 addresses cost

Global memory is accessed in 32-byte sectors (an L1/L2 cache line is 128 bytes = 4 sectors). When a warp executes a load or store, the hardware looks at the 32 lanes' addresses together and issues one transaction per distinct sector touched. That's coalescing:

• 32 lanes reading consecutive floats: 128 bytes → 4 sectors. Ideal.
• the same with a stride of 2 floats: 256 bytes spanned → 8 sectors, half of every sector wasted.
• stride ≥ 8 floats (32 bytes): every lane in its own sector → 32.
• all lanes reading the same address: 1 sector (a broadcast).

So the most important line in a kernel is the one that maps threadIdx.x to an address: consecutive lanes should touch consecutive addresses. For a row-major image that means threadIdx.x picks the column. Swap x and y and every result is still right — and every access costs 8× more.

#include <cstdint>
#include <print>
#include <set>

// Distinct 32-byte sectors touched by one warp: lane l accesses
// `bytes` bytes at byte offset addr(l).
template <class F>
int sectors(F addr, int bytes) {
    std::set<std::uint64_t> s;
    for (int lane = 0; lane < 32; ++lane)
        for (std::uint64_t b = addr(lane) / 32; b <= (addr(lane) + bytes - 1) / 32; ++b) s.insert(b);
    return static_cast<int>(s.size());
}

int main() {
    for (int stride : {1, 2, 4, 8, 32})
        std::println("float, stride {:2}: {:2} sectors", stride,
                     sectors([=](int l) { return std::uint64_t(l) * stride * 4; }, 4));
    std::println("float, offset by 1 element: {} sectors",
                 sectors([](int l) { return std::uint64_t(l + 1) * 4; }, 4));
    std::println("broadcast (same address): {} sector", sectors([](int) { return 64ull; }, 4));
    const int width = 1024;   // a row-major image of float4 pixels (16 bytes)
    std::println("float4, threadIdx.x -> column: {} sectors",
                 sectors([](int l) { return std::uint64_t(l) * 16; }, 16));
    std::println("float4, threadIdx.x -> row:    {} sectors",
                 sectors([=](int l) { return std::uint64_t(l) * width * 16; }, 16));
}

Alignment matters at the edges: cudaMalloc returns 256-byte-aligned pointers, so element 0 starts a sector; an access offset by one element touches one extra sector per warp. Structs of arrays (x[], y[], z[]) coalesce where arrays of structs often don't — the same data-oriented design argument as on the CPU, with a much bigger penalty.

5. Shared memory: tiles and bank conflicts

When a kernel needs a strided access pattern — a transpose, a column of a matrix, a stencil's neighbours — read a tile of global memory into shared memory with coalesced loads, wait, then read the tile in whatever order you need. The canonical example is the transpose:

constexpr int TILE = 32, ROWS = 8;            // block = 32 x 8 threads
__global__ void transpose(float* out, const float* in, int rows, int cols) {
    __shared__ float tile[TILE][TILE + 1];   // +1: see below
    int x = blockIdx.x * TILE + threadIdx.x;
    int y = blockIdx.y * TILE + threadIdx.y;
    for (int j = 0; j < TILE; j += ROWS)
        if (x < cols && y + j < rows)
            tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * cols + x];   // coalesced read
    __syncthreads();                          // the whole tile is in
    x = blockIdx.y * TILE + threadIdx.x;      // swap the *block* indices
    y = blockIdx.x * TILE + threadIdx.y;
    for (int j = 0; j < TILE; j += ROWS)
        if (x < rows && y + j < cols)
            out[(y + j) * rows + x] = tile[threadIdx.x][threadIdx.y + j];  // coalesced write
}

Shared memory is split into 32 banks; successive 4-byte words go to successive banks (word w is in bank w % 32). A warp's shared-memory access completes in one pass if its lanes hit 32 different banks (or the same word — that's a broadcast). Lanes hitting different words in the same bank are serialised: an n-way bank conflict costs n passes.

Reading tile[threadIdx.x][c] walks down a column. With float tile[32][32] lane l reads word 32·l + c: all 32 lanes in bank c — a 32-way conflict. Padding each row to 33 words moves lane l to bank (l + c) % 32: all different.

#include <algorithm>
#include <array>
#include <print>
#include <set>

// Passes a warp needs when lane l reads 4-byte word word(l) of shared memory.
template <class F>
int passes(F word) {
    std::array<std::set<int>, 32> per_bank;
    for (int lane = 0; lane < 32; ++lane) per_bank[word(lane) % 32].insert(word(lane));
    std::size_t worst = 0;
    for (auto& b : per_bank) worst = std::max(worst, b.size());
    return static_cast<int>(worst);
}

int main() {
    for (int pitch : {32, 33, 64, 31}) {
        std::println("pitch {:2}: row read {:2} pass(es), column read {:2} pass(es)", pitch,
                     passes([=](int l) { return 5 * pitch + l; }),   // tile[5][l]
                     passes([=](int l) { return l * pitch + 5; }));  // tile[l][5]
    }
    std::println("all lanes read word 7: {} pass (broadcast)", passes([](int) { return 7; }));
}

The rule is stated for 32-bit words; 64-bit and 128-bit accesses are split into several passes by the hardware, and Nsight Compute's shared-memory metrics tell you whether a real kernel conflicts.

Shared memory is a limited resource: 48 KB per block statically, more on recent GPUs through dynamic shared memory (extern __shared__ float buf[];, the size as the third launch parameter, and an opt-in with cudaFuncSetAttribute(..., cudaFuncAttributeMaxDynamicSharedMemorySize, bytes)). Every KB a block uses is a KB other blocks on the SM can't — see occupancy.

6. __syncthreads(): the rules

__syncthreads() is a barrier for the threads of one block: no thread continues until all have arrived, and shared- (and global-) memory writes made before it are visible to the block after it. Two rules:

1. Every thread of the block must reach the same barrier. The CUDA Programming Guide allows it in conditional code only if the condition is the same for the whole block; otherwise the kernel may hang or misbehave. if (tid < s) { ...; __syncthreads(); } is the classic bug. Let the if select who works and put the barrier after it.
2. A barrier sits between a write and another thread's read of the same shared memory (and between a read and another thread's later overwrite). Missing ones are races; extra ones only cost time.

compute-sanitizer --tool racecheck finds shared-memory races, --tool synccheck barrier misuse, --tool memcheck out-of-bounds and misaligned accesses. Run them on every kernel you change.

The block reduction is the standard exercise in both rules. Each thread puts one value in shared memory; then, for s = blockDim/2, …, 1, threads tid < s add sh[tid + s] into sh[tid], with a barrier after every step:

__global__ void block_sum(const int* in, long long* partial, int n) {
    __shared__ long long sh[256];
    int tid = threadIdx.x;
    long long acc = 0;
    for (int i = blockIdx.x * blockDim.x + tid; i < n; i += gridDim.x * blockDim.x) acc += in[i];
    sh[tid] = acc;
    __syncthreads();
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {   // same trip count for every thread
        if (tid < s) sh[tid] += sh[tid + s];
        __syncthreads();                             // outside the if
    }
    if (tid == 0) partial[blockIdx.x] = sh[0];
}

Sequential addressing (tid and tid + s) keeps the working threads contiguous, so whole warps go idle together instead of every warp diverging, and it is bank-conflict-free. Here is the same schedule with real threads and a real barrier — std::barrier plays __syncthreads():

#include <barrier>
#include <print>
#include <thread>
#include <vector>

int main() {
    constexpr int n = 8;                        // one "block" of 8 threads
    int sh[n] = {3, 1, 4, 1, 5, 9, 2, 6};       // "shared memory"
    std::barrier sync(n);
    std::vector<int> after_step;                // written by thread 0 only
    {
        std::vector<std::jthread> block;
        for (int tid = 0; tid < n; ++tid)
            block.emplace_back([&, tid] {
                for (int s = n / 2; s > 0; s /= 2) {
                    if (tid < s) sh[tid] += sh[tid + s];
                    sync.arrive_and_wait();     // every thread, every step
                    if (tid == 0) after_step.push_back(sh[0]);
                    sync.arrive_and_wait();     // nobody writes sh before 0 has read it
                }
            });
    }
    for (int v : after_step) std::print("{} ", v);
    std::println("-> sum {}", sh[0]);
}

(The second barrier in the CPU version only protects after_step's bookkeeping read; the reduction itself needs one per step.)

Modern CUDA also has warp-level primitives (__shfl_down_sync reduces within a warp without shared memory or block barriers) and cooperative groups (cg::this_thread_block().sync(), cg::reduce); learn the shared memory version first — interviews do.

7. Atomics

atomicAdd(&x, v) (and atomicMin, atomicMax, atomicCAS, …) performs an indivisible read-modify-write on global or shared memory. Use them for histograms, counters and the final combine of per-block results — not as a replacement for a reduction:

__global__ void sum_atomic(const float* in, float* total, int n) {
    // ... block reduction as above into sh[0] ...
    if (threadIdx.x == 0) atomicAdd(total, sh[0]);   // one atomic per block, not per element
}

Thousands of threads hitting one address serialise on it. And floating point atomics complete in an unspecified order, so a float sum computed with atomics can differ in the last bits from run to run — the next skill (gpu-correctness-and-profiling) is about comparing such results honestly.

8. Host ↔ device: transfers, streams, synchronisation, errors

A minimal host program, with the error checking you should never skip:

#define CUDA_CHECK(expr)                                                        \
    do {                                                                        \
        cudaError_t err_ = (expr);                                              \
        if (err_ != cudaSuccess) {                                              \
            std::fprintf(stderr, "%s:%d: %s: %s\n", __FILE__, __LINE__, #expr,  \
                         cudaGetErrorString(err_));                             \
            std::abort();                                                       \
        }                                                                       \
    } while (0)

float* d = nullptr;
CUDA_CHECK(cudaMalloc(&d, n * sizeof(float)));
CUDA_CHECK(cudaMemcpy(d, host.data(), n * sizeof(float), cudaMemcpyHostToDevice));
scale<<<(n + 255) / 256, 256>>>(d, n, 2.0f);
CUDA_CHECK(cudaGetLastError());          // launch errors: bad config, too many threads, ...
CUDA_CHECK(cudaDeviceSynchronize());     // errors *during* the kernel: illegal address, ...
CUDA_CHECK(cudaMemcpy(host.data(), d, n * sizeof(float), cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaFree(d));

• A kernel launch is asynchronous: it returns before the kernel runs and has no return value. Configuration errors (a block of 2048 threads, too much shared memory) are reported by cudaGetLastError() right after the launch. Errors inside the kernel (an out-of-bounds write that hits unmapped memory) appear at the next synchronising call — and they are sticky: the context is unusable and the process has to exit. That's why "add a cudaDeviceSynchronize() and check it" is the first step of debugging any CUDA crash.
• cudaMemcpy is synchronous with the host. cudaMemcpyAsync on a stream returns at once; it only overlaps with other work when the host buffer is pinned (cudaMallocHost / cudaHostAlloc) — from pageable memory the driver stages the copy.
• Work in one stream runs in order; work in different streams may overlap (copy engines and compute run concurrently). Synchronise with cudaStreamSynchronize(s), events (cudaEventRecord, cudaStreamWaitEvent) or cudaDeviceSynchronize(). The legacy default stream synchronises with other blocking streams; create streams with cudaStreamNonBlocking (or compile with --default-stream per-thread) when you want real concurrency.

9. Occupancy basics

An SM keeps many warps resident and switches between them every cycle to hide memory latency. Occupancy is resident warps ÷ the SM's maximum. Three per-SM budgets limit how many blocks fit: threads (warps), registers, and shared memory, plus a cap on blocks. The scarcest one wins. With the limits of compute capability 8.6 (RTX 30-series) from the Programming Guide's table — 1536 threads, 16 blocks, 64 K registers and 100 KB of shared memory per SM — and ignoring allocation granularity:

#include <algorithm>
#include <print>

int main() {
    const int sm_threads = 1536, sm_blocks = 16, sm_regs = 65536, sm_shared_kb = 100;
    struct Kernel { const char* name; int threads, regs, shared_kb; };
    for (Kernel k : {Kernel{"light", 256, 32, 0}, Kernel{"register-heavy", 256, 64, 16},
                     Kernel{"shared-heavy", 128, 32, 48}}) {
        int by_threads = sm_threads / k.threads;
        int by_regs = sm_regs / (k.regs * k.threads);
        int by_shared = k.shared_kb ? sm_shared_kb / k.shared_kb : sm_blocks;
        int blocks = std::min({by_threads, by_regs, by_shared, sm_blocks});
        std::println("{:15} blocks/SM {} (threads {}, regs {}, shared {}) -> occupancy {}%",
                     k.name, blocks, by_threads, by_regs, by_shared,
                     100 * blocks * k.threads / sm_threads);
    }
}

Higher occupancy hides latency better, up to a point; a kernel with lots of independent work per thread can be fast at 30 %. Don't guess: the CUDA runtime has cudaOccupancyMaxActiveBlocksPerMultiprocessor and cudaOccupancyMaxPotentialBlockSize, nvcc -Xptxas -v prints registers and shared memory per kernel, __launch_bounds__(256) tells the compiler the block size to budget registers for, and Nsight Compute reports achieved occupancy. Block sizes are almost always a multiple of 32 (a partial warp wastes lanes); 128–256 is the usual starting point.

10. Building CUDA code

cmake_minimum_required(VERSION 3.28)
project(splat LANGUAGES CXX CUDA)
set(CMAKE_CUDA_ARCHITECTURES native)          # or 86, "75;86;89" for a fat binary
add_library(splat_gpu STATIC composite.cu)
target_compile_features(splat_gpu PUBLIC cuda_std_20)
target_compile_options(splat_gpu PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:-lineinfo>)
find_package(CUDAToolkit REQUIRED)
target_link_libraries(splat_gpu PUBLIC CUDA::cudart)

-lineinfo maps SASS back to source lines for Nsight Compute and compute-sanitizer without slowing the code; -G builds device code unoptimised for cuda-gdb (much slower). native (CMake ≥ 3.24) compiles for the GPUs in the build machine — wrong for a CI machine without one, so list architectures explicitly there. CUDA 13 no longer targets Maxwell, Pascal or Volta GPUs; check your target hardware's compute capability before picking a toolkit. A build with no GPU present still compiles; running needs a driver — check cudaGetDeviceCount and report "no device" instead of crashing.

11. The simulator in the build task

The build task's gpu::launch runs your kernel on the CPU: blocks one after another, each block's threads as real threads that take turns from barrier to barrier, so syncthreads() really waits and every run interleaves the threads the same way. Global memory goes through gpu::global_span and shared memory through ctx.shared<T>(), and every access is recorded:

• per warp, the k-th global load (store) of every lane forms one request, costing one transaction per distinct 32-byte segment — the coalescing rule from §4;
• shared accesses are grouped the same way and cost one wavefront per access to the busiest bank — §5;
• a thread that skips a barrier its block waits at makes the block divergent (instead of hanging) — §6 rule 1; a shared byte written by one thread and touched by another between two barriers is a race — rule 2;
• out-of-range indices are refused and counted.

It's a model — no caches, no timing — but it's the same model Nsight Compute's "sectors per request" and "bank conflicts" metrics measure, and every count is deterministic, so the tests can grade them exactly.

In the project

Ungraded practice in projects/scene-viewer/gpu/ (stage 5); the evidence is your own build and test output.

1. Read gpu/preprocess_cuda.cu: one thread per Gaussian, the bounds check for the partial last block, RAII for device memory and events, cudaGetLastError after the launch. Build it with cmake --preset cuda.
2. Port the tile compositing (the per-pixel loop in core/src/render.cpp) to a kernel: one block per 16×16 tile, threadIdx.x/y → pixel, the tile's sorted splat list staged through shared memory with __syncthreads.
3. Without a GPU, the test gpu.matches_cpu is reported as skipped, not passed. With one, run it under compute-sanitizer --tool memcheck and --tool racecheck.