CUDA Programming Fundamentals

You will rarely write production CUDA. You should write enough to read it fluently, to understand what your engine is doing, and to write a custom fused op when nothing else fits.

1. What is it?#

Writing functions that run on the GPU, in C++ with CUDA extensions (or in Triton, which is easier and covered in Section XIII.11).

2. Why learn it?#

Three concrete reasons:

To read kernels. vLLM’s csrc/, FlashAttention’s implementation, CUTLASS — the source of truth for how things actually work.
To write fused ops when your architecture has something the library doesn’t cover.
To interpret profiles. Nsight’s output is meaningless without knowing what a kernel does.

3. Your first kernel#

// vecadd.cu
#include <cstdio>

__global__ void vecadd(const float* a, const float* b, float* c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) c[i] = a[i] + b[i];
}

int main() {
    int n = 1 << 20;
    size_t bytes = n * sizeof(float);

    float *h_a = (float*)malloc(bytes), *h_b = (float*)malloc(bytes),
          *h_c = (float*)malloc(bytes);
    for (int i = 0; i < n; i++) { h_a[i] = i; h_b[i] = 2*i; }

    float *d_a, *d_b, *d_c;
    cudaMalloc(&d_a, bytes); cudaMalloc(&d_b, bytes); cudaMalloc(&d_c, bytes);
    cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);

    int threads = 256, blocks = (n + threads - 1) / threads;
    vecadd<<<blocks, threads>>>(d_a, d_b, d_c, n);
    cudaDeviceSynchronize();

    cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
    printf("c[42] = %f (expect %f)\n", h_c[42], 3.0f*42);

    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
    free(h_a); free(h_b); free(h_c);
}

nvcc -O3 -arch=sm_80 vecadd.cu -o vecadd && ./vecadd

Always check errors. The single most common CUDA beginner mistake is ignoring return codes, because errors are asynchronous and surface much later at a confusing place:

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

CUDA_CHECK(cudaMalloc(&d_a, bytes));
vecadd<<<blocks, threads>>>(d_a, d_b, d_c, n);
CUDA_CHECK(cudaGetLastError());        // catches launch errors
CUDA_CHECK(cudaDeviceSynchronize());   // catches execution errors

4. The essential language features#

__global__   runs on GPU, called from host. Must return void.
__device__   runs on GPU, called from GPU.
__host__     runs on CPU (default).
__shared__   shared memory, block-scoped
__constant__ read-only, broadcast-optimized (64 KB)
__restrict__ promises no aliasing — enables optimization

Built-in variables:
  threadIdx.{x,y,z}   thread index within the block
  blockIdx.{x,y,z}    block index within the grid
  blockDim.{x,y,z}    block size
  gridDim.{x,y,z}     grid size

Synchronization:
  __syncthreads()             barrier for all threads in the block
  __syncwarp(mask)            barrier within a warp
  __threadfence()             memory ordering

Warp primitives (very fast — data is already co-located):
  __shfl_sync(mask, val, srcLane)
  __shfl_down_sync(mask, val, delta)
  __ballot_sync(mask, pred)
  __any_sync / __all_sync
  __reduce_add_sync (Ampere+)

Atomics:
  atomicAdd, atomicMax, atomicCAS, ...

5. Kernel 2: a reduction (the fundamental pattern)#

Summing an array teaches warp primitives, shared memory, and the two-level reduction pattern that appears in every norm, softmax, and attention kernel.

__inline__ __device__ float warpReduceSum(float val) {
    for (int offset = 16; offset > 0; offset /= 2)
        val += __shfl_down_sync(0xffffffff, val, offset);
    return val;                       // lane 0 has the warp's sum
}

__inline__ __device__ float blockReduceSum(float val) {
    static __shared__ float shared[32];      // one slot per warp
    int lane = threadIdx.x % 32;
    int wid  = threadIdx.x / 32;

    val = warpReduceSum(val);                 // reduce within each warp
    if (lane == 0) shared[wid] = val;         // one value per warp to shared
    __syncthreads();

    // first warp reduces the per-warp partial sums
    val = (threadIdx.x < blockDim.x / 32) ? shared[lane] : 0.0f;
    if (wid == 0) val = warpReduceSum(val);
    return val;                               // thread 0 has the block's sum
}

__global__ void reduce(const float* in, float* out, int n) {
    float sum = 0.0f;
    // grid-stride loop: handles any n with any grid size
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n;
         i += blockDim.x * gridDim.x)
        sum += in[i];
    sum = blockReduceSum(sum);
    if (threadIdx.x == 0) atomicAdd(out, sum);
}

Three patterns to internalize here:

Grid-stride loop — decouples grid size from problem size; the same kernel works for any n.
Warp shuffle reduction — no shared memory, no barriers, ~5 instructions.
Two-level reduction — warp then block. Every norm and softmax kernel does this.

6. Kernel 3: tiled matmul#

The canonical shared-memory kernel.

#define TILE 32

__global__ void matmul(const float* A, const float* B, float* C, int M, int N, int K) {
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE + 1];      // +1 to avoid bank conflicts

    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float acc = 0.0f;

    for (int t = 0; t < (K + TILE - 1) / TILE; t++) {
        // cooperatively load one tile of A and one of B
        int aCol = t * TILE + threadIdx.x;
        int bRow = t * TILE + threadIdx.y;
        As[threadIdx.y][threadIdx.x] = (row < M && aCol < K) ? A[row*K + aCol] : 0.0f;
        Bs[threadIdx.y][threadIdx.x] = (bRow < K && col < N) ? B[bRow*N + col] : 0.0f;
        __syncthreads();

        #pragma unroll
        for (int k = 0; k < TILE; k++)
            acc += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        __syncthreads();                       // before overwriting the tiles
    }
    if (row < M && col < N) C[row*N + col] = acc;
}

This achieves maybe 20-30% of cuBLAS. The gap comes from: no tensor cores, no double buffering, one output element per thread (should be 4-8), no vectorized loads, no async copy. Closing that gap is a multi-week project — which is why you use cuBLAS.

But writing this once teaches you what every GEMM kernel is doing, which is the point.

7. Integrating with PyTorch#

// my_op.cu
#include <torch/extension.h>

__global__ void scale_kernel(float* x, float s, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) x[i] *= s;
}

void scale(torch::Tensor x, float s) {
    TORCH_CHECK(x.is_cuda(), "x must be on CUDA");
    TORCH_CHECK(x.is_contiguous(), "x must be contiguous");
    int n = x.numel();
    int threads = 256, blocks = (n + threads - 1) / threads;
    scale_kernel<<<blocks, threads>>>(x.data_ptr<float>(), s, n);
    C10_CUDA_KERNEL_LAUNCH_CHECK();
}

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { m.def("scale", &scale); }

from torch.utils.cpp_extension import load
mod = load(name="my_op", sources=["my_op.cu"], verbose=True)
x = torch.ones(1000, device='cuda')
mod.scale(x, 3.0)

torch.utils.cpp_extension.load compiles on the fly — good for development. For production, use setup.py with CUDAExtension and ship a built wheel.

8. Debugging#

# Memory errors (out of bounds, misaligned, leaks)
compute-sanitizer --tool memcheck ./myapp

# Race conditions on shared memory
compute-sanitizer --tool racecheck ./myapp

# Uninitialized memory
compute-sanitizer --tool initcheck ./myapp

# Synchronization errors
compute-sanitizer --tool synccheck ./myapp

# printf works in kernels (slow, but works)
printf("thread %d: val=%f\n", threadIdx.x, val);

# Debug build with device-side debugging
nvcc -g -G ...     # then cuda-gdb

compute-sanitizer is your best friend. Run it on any new kernel. Out-of-bounds writes on a GPU corrupt memory silently and manifest as wrong answers somewhere else entirely.

9. Common mistakes#

No bounds check. if (i < n) — always.

Forgetting __syncthreads() after writing shared memory. Intermittent wrong answers.

__syncthreads() inside divergent code. Undefined behavior.

Ignoring error codes. Errors surface later, confusingly.

Assuming threads in a block execute simultaneously. They’re warps, scheduled independently.

Bank conflicts from unpadded [32][32] shared arrays.

Integer overflow in index arithmetic. blockIdx.x * blockDim.x can overflow int for large grids; use size_t or int64_t.

Race on atomicAdd initialization. Zero the output before launching.

10. Hands-on exercise#

A. Vector add. Write, compile, run, and verify. Sweep block size and plot achieved bandwidth. What fraction of peak do you get? This is Project 03.

B. Reduction. Implement the two-level reduction. Verify against torch.sum. Benchmark against it. How close do you get?

C. Tiled matmul. Implement it. Verify correctness. Benchmark against cuBLAS for N=1024, 2048, 4096. Compute your achieved TFLOP/s and the fraction of peak.

D. Fused RMSNorm. Write a fused RMSNorm kernel (reduction + scale) and integrate it with PyTorch. Benchmark against the unfused PyTorch version. This is a genuinely useful kernel.

E. Sanitize. Deliberately introduce an out-of-bounds write and confirm compute-sanitizer catches it. Then introduce a shared-memory race and catch it with racecheck.

F. Read real code. Read vllm/csrc/layernorm_kernels.cu and activation_kernels.cu. They’re short. Identify: the grid-stride loop, the reduction, the vectorized loads.

11. Interview questions#

Write a vector-add kernel from memory, including the bounds check.
Explain the two-level reduction pattern and why warp shuffles are used.
What is a grid-stride loop and why use it?
Why does __shared__ float tile[32][33] have a 33?
Why must all threads in a block reach __syncthreads()?
How would you debug a kernel producing wrong answers intermittently?
Why does a hand-written tiled matmul only reach 30% of cuBLAS?

12. Further reading#

[FUNDAMENTAL] Kirk & Hwu, Programming Massively Parallel Processors (4th ed.)
[REFERENCE] CUDA C++ Programming Guide and Best Practices Guide
[REFERENCE] NVIDIA cuda-samples repository
[ESTABLISHED] “How to Optimize a CUDA Matmul Kernel” (siboehm) — an excellent step-by-step walkthrough from naive to near-cuBLAS
Next: 09 — Coalescing and access patterns