Learn: Your First CUDA Kernels

What you're learning

How to actually write the code that runs on the hardware from Lessons 01–03. By the end you’ll have written, compiled, and timed two real kernels — vector add and a naive matmul — and you’ll be able to explain, from the roofline, why the naive matmul wastes >99% of the GPU. That “why” is the entire motivation for Lesson 05.

This is Lesson 04 of the GPU-systems curriculum. It assumes Tier 1: 01-gpu-architecture-and-simt (SIMT, warps, the grid/block/thread mapping), 02-gpu-memory-hierarchy (registers → SMEM → HBM, coalescing), and 03-performance-modeling-roofline (arithmetic intensity, the ridge point). We now stop describing the machine and start programming it.

1. Learning map

graph TD
    A["Lessons 01–03<br/>SIMT · memory · roofline"] --> B["CUDA programming model<br/>host vs device, kernels"]
    B --> C["Qualifiers<br/>__global__ / __device__ / __host__"]
    B --> D["Launch config<br/>&lt;&lt;&lt;grid, block&gt;&gt;&gt;"]
    D --> E["Global thread index<br/>idx = blockIdx·blockDim + threadIdx"]
    E --> F["Boundary guard<br/>if (idx &lt; n)"]
    B --> G["Memory management<br/>cudaMalloc / Memcpy / Free"]
    G --> H["Unified Memory<br/>cudaMallocManaged"]
    E --> I["Vector add<br/>(end-to-end)"]
    F --> I
    G --> I
    I --> J["Grid-stride loops<br/>arbitrary n, decoupled from grid"]
    I --> K["Error checking<br/>CUDA_CHECK + getLastError"]
    K --> L["Async launch<br/>cudaDeviceSynchronize"]
    L --> M["CUDA events<br/>timing"]
    I --> N["Naive matmul<br/>one thread / output elem"]
    N --> O["WHY it's slow:<br/>redundant row/col HBM reads<br/>I ≈ 0.25 FLOP/byte"]
    O --> P["Roofline: memory-bound<br/>→ Lesson 05 tiling"]

    style B fill:#44a,color:#fff
    style I fill:#4a4,color:#fff
    style O fill:#a44,color:#fff

The spine: the programming model → how a launch turns into threads → how a thread finds its data → a full working kernel → the robustness/timing scaffolding around it → a second kernel that is correct but slow, and the roofline diagnosis of why.

2. Why this matters

You can recite the SIMT model and the roofline all day, but a GPU researcher’s day job is writing and reading kernels. Every abstraction in the previous three lessons — warps, coalescing, occupancy, arithmetic intensity — only becomes actionable once you can express work as a __global__ function and reason about the exact bytes each thread touches.

Vector add is the “hello world” that establishes the mechanical skeleton every CUDA program shares: allocate on the device, copy in, launch, copy out, free. Naive matmul is the “hello world” of why GPUs are hard: it’s trivially correct, embarrassingly parallel, maps perfectly onto SIMT — and still runs at a few percent of peak because it drowns in redundant HBM traffic. Understanding that gap mechanically is the difference between someone who writes CUDA and someone who optimizes it. This lesson gets you the first; it sets up everything that gets you the second.

The through-line to LLMs

A transformer is a chain of GEMMs (Lesson 01 §9). The naive matmul here is the pedagogical ancestor of cuBLAS/CUTLASS GEMM and FlashAttention. Every optimization in Lessons 05–06 is a scheme to fix the exact pathology you’ll diagnose in §9: too little work done per byte pulled from HBM.

3. The CUDA programming model: host vs device

An unconditional truth to anchor on: the CPU and GPU are separate processors with separate memories, connected by a bus (PCIe/NVLink). The CPU is the host; the GPU is the device. They do not share an address space by default — a raw pointer you got from malloc is a host pointer and dereferencing it on the device is undefined, and vice versa. Almost every beginner bug traces back to forgetting which side of the bus a pointer lives on.

The host orchestrates; the device computes. A typical program is:

sequenceDiagram
    participant H as Host (CPU)
    participant D as Device (GPU)
    H->>D: cudaMalloc — reserve device memory
    H->>D: cudaMemcpy H2D — copy inputs over the bus
    H->>D: kernel<<<grid, block>>>(...)  (async)
    Note over D: thousands of threads run the kernel
    D-->>H: cudaMemcpy D2H — copy results back
    H->>D: cudaFree — release device memory

A kernel is a function you write once, from the point of view of a single thread, that the hardware then runs across a whole grid of threads in parallel (the SIMT model of Lesson 01). You do not write a loop over threads — you write the body of one thread and describe how many to launch.

Function qualifiers

CUDA C++ adds three execution-space qualifiers that tell nvcc where a function runs and from where it can be called:¹

Qualifier	Runs on	Callable from	Use
`__global__`	device	host (and device, dynamic parallelism)	a kernel — the launch entry point. Must return `void`.
`__device__`	device	device only	helper functions called inside a kernel
`__host__`	host	host	ordinary CPU function (the default if unqualified)

__host__ __device__ float square(float x) { return x * x; }  // compiled for BOTH
__device__        float relu(float x) { return x > 0.f ? x : 0.f; }  // device-only helper
 
__global__ void apply(float* y, const float* x, int n) {   // a kernel
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) y[idx] = relu(square(x[idx]));            // calls __device__ helpers
}

__global__ vs __device__ is the distinction people fumble

__global__ is a kernel — it is launched with the <<<...>>> syntax and always returns void. __device__ is a helper — it is called like a normal function, but only from device code, and it can return a value. You call __device__ functions from inside __global__ (or other __device__) functions. You never launch a __device__ function and you never call a __global__ function like a normal function.

4. Launch configuration and the global thread index

You launch a kernel with the triple-angle-bracket syntax:

kernel<<<gridDim, blockDim>>>(arg1, arg2, ...);

blockDim — how many threads per block. This is a dim3 (up to 3D). A block runs entirely on one SM (Lesson 01 §6), so its size is capped at 1024 threads.¹ Common choice: 128 or 256.
gridDim — how many blocks. Also a dim3. This is how you scale past one SM to the whole GPU; grids routinely have millions of blocks.

Inside the kernel, every thread gets read-only built-in variables identifying which thread it is:

threadIdx.{x,y,z} — this thread’s index within its block.
blockIdx.{x,y,z} — this block’s index within the grid.
blockDim.{x,y,z} — the block’s dimensions (same for all threads).
gridDim.{x,y,z} — the grid’s dimensions.

The fundamental job of the first lines of nearly every kernel is to turn these local coordinates into a single global index that says which piece of data this thread owns. For a 1D launch:

int idx = blockIdx.x * blockDim.x + threadIdx.x;

Read it as: skip all the threads in the blocks before me (blockIdx.x * blockDim.x), then add my offset within my own block (threadIdx.x). This gives every thread in the entire grid a unique, contiguous index — and, crucially, consecutive threads get consecutive indices, which (Lesson 02 §6) is exactly what makes global memory accesses coalesce.

The thread-index → data-element mapping (vector add)

idx = blockIdx.x · blockDim.x + threadIdx.x

Array C[idx] = A[idx] + B[idx]

0123
4567
891011

Consecutive threads → consecutive indices → coalesced 128-byte loads from HBM.

Figure 1 — Each thread computes exactly one output element. With blockDim.x = 4 and 3 blocks, thread (blockIdx=1, threadIdx=2) maps to idx = 1·4 + 2 = 6. Because lane order is preserved, a warp’s 32 threads read 32 contiguous words — one coalesced transaction (Lesson 02 §6).

The boundary guard

You will almost never have n be an exact multiple of blockDim. You compute the number of blocks by rounding up:

int threadsPerBlock = 256;
int blocks = (n + threadsPerBlock - 1) / threadsPerBlock;   // ceil(n / 256)

This launches at least n threads — but usually a few more. The extra threads in the last block have idx >= n and must do nothing, or they’ll read/write out of bounds. Hence the ubiquitous guard:

if (idx < n) { /* real work */ }

Omitting the boundary guard is a silent memory-corruption bug

Without if (idx < n), the surplus threads in the final block index past the end of your arrays. On the write side that’s an out-of-bounds store — it can clobber unrelated device memory or trip cudaErrorIllegalAddress (which, being async, may only surface on a later API call, making it maddening to localize). On the read side you get garbage. It often “works” in tiny tests where n happens to divide evenly, then explodes in production. Always guard, and always run under compute-sanitizer (§8).

Quiz: Launch config & the boundary guard. You have n = 1000 elements and pick threadsPerBlock = 256. (a) How many blocks does (n + tpb - 1)/tpb give, and how many threads launch in total? (b) How many threads have idx >= n, and what do they do? (c) What breaks if you instead compute blocks = n / tpb?

Answer

(a) $(1000 + 255) /256 = 1255/256 = 4$ blocks (integer division), $4 \times 256 = 1024$ threads launched. The ceiling-division idiom guarantees at least n threads.
(b) $1024 - 1000 = 24$ threads have idx in $[1000, 1023]$ , i.e. idx >= n. The if (idx < n) guard makes them do nothing — no load, no store. They still occupy warp slots (the last block’s last warp is partly idle), a tiny, unavoidable overhead.
(c) blocks = 1000/256 = 3 (truncating) → only $3 \times 256 = 768$ threads → elements 768..999 are never computed. Their C entries keep whatever garbage was in device memory. Silent wrong answers on the tail — worse than a crash, because it looks like it worked. Always round up.

5. Memory management

Because host and device memories are separate, you explicitly manage the device side. Four core calls:

float *d_x = nullptr;
size_t bytes = n * sizeof(float);
 
cudaMalloc(&d_x, bytes);                         // allocate on device (note: &d_x)
cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice);   // H2D: copy input in
// ... launch kernel that reads/writes d_x ...
cudaMemcpy(h_x, d_x, bytes, cudaMemcpyDeviceToHost);   // D2H: copy result out
cudaFree(d_x);                                   // release

Points that trip people up:

cudaMalloc(&ptr, bytes) takes the address of your pointer. It writes the new device address into your pointer variable, so it needs void**. Passing ptr instead of &ptr is a classic bug.
The direction flag must match reality. cudaMemcpyHostToDevice (H2D) copies host→device; cudaMemcpyDeviceToHost (D2H) the reverse. There’s also DeviceToDevice, and cudaMemcpyDefault which infers direction from the pointers (works with Unified Virtual Addressing).
cudaMemcpy is synchronous (for pageable host memory): it blocks the host until the copy completes, and implicitly waits for prior device work on the default stream. That’s actually convenient here — it means the D2H copy after a kernel launch won’t run until the kernel finishes (more on this in §8).
A device pointer is not dereferenceable on the host. d_x[0] in host code is undefined behavior. You must copy back first.

Unified Memory: `cudaMallocManaged`

There’s a friendlier alternative that removes the explicit copies:

float *x;
cudaMallocManaged(&x, bytes);   // one pointer, valid on BOTH host and device
for (int i = 0; i < n; i++) x[i] = i;   // touch on host
kernel<<<blocks, tpb>>>(x, n);          // touch on device — no cudaMemcpy needed
cudaDeviceSynchronize();                // must sync before reading on host again
printf("%f\n", x[0]);                   // read on host
cudaFree(x);

Unified Memory (UM) gives you a single pointer whose backing pages the driver migrates on demand between host and device (page-faulting them across the bus as each side touches them).¹ It’s great for prototyping, complex pointer-based data structures, and oversubscription (allocating more than device memory). The cost: those migrations are implicit HBM/PCIe traffic you don’t see in your code, so for performance-critical paths you usually want explicit cudaMalloc + cudaMemcpy (and cudaMemPrefetchAsync / cudaMemAdvise when using UM seriously) so the data movement is under your control. We use explicit memory below because being deliberate about bytes crossing the bus is the discipline this curriculum is about.

Managed memory is a convenience, not a free lunch

cudaMallocManaged doesn’t make data movement disappear — it hides it behind page faults. On the roofline (Lesson 03), those faults are still HBM/PCIe bytes; you’ve just lost visibility into them. Reach for it to get something working fast, then switch to explicit copies (or prefetch hints) when you start counting bytes.

6. Vector add — end to end

Here is the complete, compilable program. Save as vadd.cu.

// vadd.cu — the canonical first CUDA program: C = A + B
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
 
// ---- error-checking macro (explained in §8) ----
#define CUDA_CHECK(call)                                                      \
    do {                                                                      \
        cudaError_t err_ = (call);                                            \
        if (err_ != cudaSuccess) {                                            \
            fprintf(stderr, "CUDA error %s:%d: '%s' -> %s\n",                 \
                    __FILE__, __LINE__, #call, cudaGetErrorString(err_));     \
            exit(EXIT_FAILURE);                                               \
        }                                                                     \
    } while (0)
 
// ---- the kernel: one thread per element ----
__global__ void vecAdd(const float* A, const float* B, float* C, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {                 // boundary guard — see §4
        C[idx] = A[idx] + B[idx];
    }
}
 
int main() {
    const int n = 1 << 20;                 // ~1M elements (not a multiple-only case)
    const size_t bytes = n * sizeof(float);
 
    // 1. allocate + init on host
    float *h_A = (float*)malloc(bytes);
    float *h_B = (float*)malloc(bytes);
    float *h_C = (float*)malloc(bytes);
    for (int i = 0; i < n; i++) { h_A[i] = 1.0f; h_B[i] = 2.0f; }
 
    // 2. allocate on device
    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc(&d_A, bytes));
    CUDA_CHECK(cudaMalloc(&d_B, bytes));
    CUDA_CHECK(cudaMalloc(&d_C, bytes));
 
    // 3. copy inputs H2D
    CUDA_CHECK(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice));
 
    // 4. launch
    int threadsPerBlock = 256;
    int blocks = (n + threadsPerBlock - 1) / threadsPerBlock;   // ceil div
    vecAdd<<<blocks, threadsPerBlock>>>(d_A, d_B, d_C, n);
    CUDA_CHECK(cudaGetLastError());        // catches launch-config errors (§8)
    CUDA_CHECK(cudaDeviceSynchronize());   // wait + catch async runtime errors
 
    // 5. copy result D2H
    CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));
 
    // 6. verify
    double maxErr = 0.0;
    for (int i = 0; i < n; i++) maxErr = fmax(maxErr, fabs(h_C[i] - 3.0f));
    printf("n = %d, max error = %f\n", n, maxErr);   // expect 0.000000
 
    // 7. cleanup
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    free(h_A); free(h_B); free(h_C);
    return 0;
}

Compile and run:

nvcc -O3 -arch=sm_90 vadd.cu -o vadd   # sm_90 = Hopper/H100; use sm_80 for A100
./vadd
# -> n = 1048576, max error = 0.000000

Every CUDA program you ever write is a variation on this skeleton: allocate → copy in → launch → sync → copy out → free. Memorize its shape.²

Vector add is not a good benchmark of the GPU

This kernel does 1 FLOP per 3 elements touched (read A, read B, write C = 12 bytes), so its arithmetic intensity is $I = 1/12 \approx 0.08$ FLOP/byte — off the far-left, deeply memory-bound end of the roofline (Lesson 03). It runs at HBM bandwidth and touches ~0% of the FP32 ALUs, let alone the Tensor Cores. It’s the right first program (it teaches the mechanics) and the wrong thing to be impressed by (it tells you nothing about compute throughput). That’s precisely why §9’s matmul matters.

7. Grid-stride loops

The vector-add kernel above assumes you launch at least n threads. That couples your grid size to your data size, and it breaks or wastes launches when n is huge or unknown. The idiomatic fix is the grid-stride loop: launch a fixed, hardware-sized grid, and have each thread walk over multiple elements, striding by the total number of threads in the grid.³

__global__ void vecAddGridStride(const float* A, const float* B, float* C, int n) {
    int stride = gridDim.x * blockDim.x;                     // total threads in grid
    for (int idx = blockIdx.x * blockDim.x + threadIdx.x;    // this thread's start
         idx < n;
         idx += stride) {                                    // jump a whole grid
        C[idx] = A[idx] + B[idx];
    }
}

If the grid has T total threads and n > T, thread k handles elements k, k+T, k+2T, …. The idx < n condition doubles as both the loop bound and the boundary guard — one construct handles arbitrary n.

Why prefer this over “just launch n threads”?

Decouples grid size from problem size. You can launch a grid sized to the GPU (e.g. blocks = 32 * numSMs) and it processes any n, including n larger than the max grid dimension (~2³¹−1 blocks in x, but you rarely want that many tiny blocks).¹
Coalescing is preserved. Because the stride is the whole grid width, each iteration a warp still reads 32 contiguous elements — every pass is coalesced (Lesson 02 §6). Contrast with the naive alternative of giving each thread a contiguous chunk (idx*K … idx*K+K), which strides within a warp and destroys coalescing.
Amortizes per-thread setup and enables reuse. The index math and any per-thread constants are computed once, then reused across many elements — useful when the body is more than a single add.
Tunable occupancy & debuggability. You can launch a single block of a single thread and the loop still processes the entire array correctly — invaluable for debugging. Then scale the grid up for performance without touching kernel logic.

The launch just picks a grid sized to the device rather than to n:

int tpb = 256;
int numSMs;
cudaDeviceGetAttribute(&numSMs, cudaDevAttrMultiProcessorCount, 0);
int blocks = 32 * numSMs;              // enough blocks to fill the GPU with headroom
vecAddGridStride<<<blocks, tpb>>>(d_A, d_B, d_C, n);   // works for ANY n

Quiz: Grid-stride loops. Why does striding by gridDim.x * blockDim.x keep memory accesses coalesced, whereas giving each thread a contiguous block of K elements (for j in [idx*K, idx*K+K)) does not?

Answer

Coalescing is a per-warp property (Lesson 02 §6): the 32 lanes of a warp issue one memory instruction, and it’s fast only if their 32 addresses fall in one/few 128-byte lines — i.e. the lanes access consecutive words.

Grid-stride: on each iteration, lane t of a warp accesses base + t for the same base (the warp moves as a unit, then all lanes jump by stride). So within any single iteration the 32 lanes touch 32 contiguous words → one coalesced 128-byte transaction. Every pass stays coalesced.

Contiguous chunks: lane t starts at t*K, lane t+1 at (t+1)*K — the lanes are K elements apart within the same instruction. For K larger than 1 that’s a strided access spanning up to 32 different cache lines → up to a 32× effective-bandwidth loss. The chunking “feels” cache-friendly from a single thread’s CPU-style view, but it’s exactly backwards for SIMT: you want neighboring threads, not neighboring iterations of one thread, to hit neighboring addresses.

8. Error checking, async launches, and synchronization

CUDA’s C API returns error codes; it does not throw. And kernel launches are asynchronous — the launch call returns to the host immediately, before the kernel has finished (often before it has even started).¹ These two facts make robust error handling non-obvious, and getting it wrong is the most common reason a beginner’s kernel “does nothing” or “sometimes works.”

The `CUDA_CHECK` macro

Wrap every runtime API call. The macro from §6, repeated for emphasis:

#define CUDA_CHECK(call)                                                      \
    do {                                                                      \
        cudaError_t err_ = (call);                                            \
        if (err_ != cudaSuccess) {                                            \
            fprintf(stderr, "CUDA error %s:%d: '%s' -> %s\n",                 \
                    __FILE__, __LINE__, #call, cudaGetErrorString(err_));     \
            exit(EXIT_FAILURE);                                               \
        }                                                                     \
    } while (0)

The do { … } while(0) makes the macro a single statement that behaves correctly after an if without braces. #call stringifies the call so the message names the exact failing line.

Checking a kernel launch — the two-step

You cannot wrap kernel<<<...>>>(...) in CUDA_CHECK directly (the launch syntax isn’t an expression returning cudaError_t). And a launch can fail in two distinct ways, caught differently:

myKernel<<<blocks, tpb>>>(args...);
CUDA_CHECK(cudaGetLastError());       // (1) SYNCHRONOUS launch errors:
                                      //     bad config (too many threads/block,
                                      //     too much shared mem), invalid args
CUDA_CHECK(cudaDeviceSynchronize());  // (2) waits for the kernel, then reports
                                      //     ASYNC runtime errors (illegal address,
                                      //     misaligned access) raised DURING execution

cudaGetLastError() (right after launch) catches errors detected at launch time — e.g. blockDim > 1024, or requesting more shared memory than the SM has. It also clears the sticky error state.
cudaDeviceSynchronize() blocks the host until all prior device work completes, and returns any error that occurred during kernel execution (like an out-of-bounds access from a missing boundary guard). Because the launch is async, such an error would otherwise surface only on some later unrelated API call, pointing you at the wrong line.

The async-launch trap: a kernel error can "teleport" to a later call

Kernel launches return before the kernel runs. If a kernel does something illegal (e.g. an out-of-bounds write from a missing if (idx < n)), the error isn’t reported by the launch itself — it’s sticky and gets returned by the next synchronizing call (a cudaMemcpy, a cudaDeviceSynchronize, or even the next launch). Without an explicit cudaGetLastError() + cudaDeviceSynchronize() right after each launch, your program will blame the wrong line — often a cudaMemcpy that is itself perfectly fine. This is why the two-step check exists.

Don't leave cudaDeviceSynchronize() in hot paths in production

The two-step check is essential during development, but cudaDeviceSynchronize() serializes the host and device (kills overlap and pipelining). In release builds, gate the sync behind a debug flag, or rely on the fact that the next real synchronizing call (a cudaMemcpy, stream sync, or event query) will surface a sticky error anyway. During bring-up: always sync-and-check. In a throughput-critical loop: check cudaGetLastError() cheaply and sync sparingly.

compute-sanitizer is your valgrind for kernels

Run compute-sanitizer ./vadd (ships with the toolkit; formerly cuda-memcheck).⁴ It catches out-of-bounds accesses, misaligned loads, race conditions on shared memory, and uninitialized reads — pinpointing the exact kernel and thread. Make it a reflex before trusting a new kernel’s output.

Quiz: Async launch semantics. A colleague times a kernel with auto t0 = clock(); kernel<<<g,b>>>(...); auto t1 = clock(); and reports "3 microseconds — blazing fast!" Separately, their program crashes in a cudaMemcpy that they've triple-checked is correct. Explain both, and fix both.

Answer

Both symptoms are the same root cause: kernel launches are asynchronous — the launch call returns to the host almost immediately, before (often long before) the kernel finishes.

The “3 µs” measurement is just the launch overhead — the time to enqueue the kernel into the stream, not to execute it. The host raced ahead to t1 while the GPU was still working. Fix: time with CUDA events recorded in the stream (cudaEventRecord around the launch + cudaEventSynchronize + cudaEventElapsedTime), which measure the actual device execution window.

The “crash in a correct cudaMemcpy” is a sticky error teleport: the kernel actually did something illegal (e.g. an out-of-bounds access), but because the launch already returned, that error wasn’t reported until the next synchronizing call — the cudaMemcpy — which returns the earlier kernel’s error code. The cudaMemcpy is innocent; it’s just the messenger. Fix: put CUDA_CHECK(cudaGetLastError()); CUDA_CHECK(cudaDeviceSynchronize()); immediately after the launch so the error is reported at its true source, then run compute-sanitizer to pinpoint the offending thread.

Timing with CUDA events

Because launches are async, you cannot time a kernel with CPU std::chrono around the launch line — you’d measure the launch overhead, not the execution. The correct tool is CUDA events, which are timestamps recorded in the device’s stream:¹

cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
 
// warm-up launch (first launch pays JIT/context init costs — don't time it)
vecAdd<<<blocks, tpb>>>(d_A, d_B, d_C, n);
CUDA_CHECK(cudaDeviceSynchronize());
 
CUDA_CHECK(cudaEventRecord(start));           // timestamp in the stream
vecAdd<<<blocks, tpb>>>(d_A, d_B, d_C, n);
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));       // wait until 'stop' has been reached
 
float ms = 0.f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));  // milliseconds, ~0.5us resolution
 
// turn it into achieved bandwidth (vector add moves 3 arrays):
double gb = 3.0 * n * sizeof(float) / 1e9;
printf("vecAdd: %.3f ms, %.1f GB/s\n", ms, gb / (ms / 1e3));

Events sit in the stream, so they bracket exactly the device work between them. Compare the printed GB/s to your GPU’s spec HBM bandwidth (H100 ≈ 3.35 TB/s)⁵ — for a well-behaved memory-bound kernel like vector add you should reach a large fraction of it, which is itself the confirmation that vector add is bandwidth-bound.

9. Naive matmul — correct, and instructively slow

Now the workload that actually matters for LLMs: $C = A \cdot B$ with $A \in R^{M \times K}$ , $B \in R^{K \times N}$ , $C \in R^{M \times N}$ . The obvious parallelization: one thread per output element $C_{ij}$ , using a 2D grid.²

// matmul_naive.cu — one thread computes one output element C[row][col]
#include <cstdio>
#include <cuda_runtime.h>
 
#define CUDA_CHECK(call)                                                      \
    do { cudaError_t e_=(call); if(e_!=cudaSuccess){                          \
        fprintf(stderr,"CUDA %s:%d %s -> %s\n",__FILE__,__LINE__,#call,       \
                cudaGetErrorString(e_)); exit(1);} } while(0)
 
// Row-major: A is M x K, B is K x N, C is M x N.  A[i][k] = A[i*K + k].
__global__ void matmulNaive(const float* A, const float* B, float* C,
                            int M, int N, int K) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;   // -> N dimension
    int row = blockIdx.y * blockDim.y + threadIdx.y;   // -> M dimension
    if (row < M && col < N) {                          // 2D boundary guard
        float acc = 0.0f;
        for (int k = 0; k < K; ++k) {
            acc += A[row * K + k] * B[k * N + col];     // full row of A · full col of B
        }
        C[row * N + col] = acc;
    }
}
 
int main() {
    int M = 1024, N = 1024, K = 1024;
    size_t bA = (size_t)M*K*sizeof(float);
    size_t bB = (size_t)K*N*sizeof(float);
    size_t bC = (size_t)M*N*sizeof(float);
 
    float *h_A=(float*)malloc(bA), *h_B=(float*)malloc(bB), *h_C=(float*)malloc(bC);
    for (int i=0;i<M*K;i++) h_A[i]=1.0f;
    for (int i=0;i<K*N;i++) h_B[i]=1.0f;               // so every C[i][j] == K
 
    float *d_A,*d_B,*d_C;
    CUDA_CHECK(cudaMalloc(&d_A,bA)); CUDA_CHECK(cudaMalloc(&d_B,bB)); CUDA_CHECK(cudaMalloc(&d_C,bC));
    CUDA_CHECK(cudaMemcpy(d_A,h_A,bA,cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B,h_B,bB,cudaMemcpyHostToDevice));
 
    dim3 block(16, 16);                                // 256 threads, 2D
    dim3 grid((N + block.x - 1)/block.x,               // ceil(N/16) in x
              (M + block.y - 1)/block.y);              // ceil(M/16) in y
 
    // time it
    cudaEvent_t s,e; CUDA_CHECK(cudaEventCreate(&s)); CUDA_CHECK(cudaEventCreate(&e));
    matmulNaive<<<grid,block>>>(d_A,d_B,d_C,M,N,K);    // warm-up
    CUDA_CHECK(cudaGetLastError()); CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaEventRecord(s));
    matmulNaive<<<grid,block>>>(d_A,d_B,d_C,M,N,K);
    CUDA_CHECK(cudaEventRecord(e)); CUDA_CHECK(cudaEventSynchronize(e));
    float ms=0; CUDA_CHECK(cudaEventElapsedTime(&ms,s,e));
 
    CUDA_CHECK(cudaMemcpy(h_C,d_C,bC,cudaMemcpyDeviceToHost));
    double gflops = 2.0*M*N*K / (ms/1e3) / 1e9;
    printf("C[0]=%.1f (expect %d)  |  %.3f ms  |  %.1f GFLOP/s\n", h_C[0], K, ms, gflops);
 
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C); free(h_A); free(h_B); free(h_C);
    return 0;
}

nvcc -O3 -arch=sm_90 matmul_naive.cu -o matmul_naive && ./matmul_naive

It’s correct — every $C_{ij}$ comes out to $K$ . It’s also embarrassingly parallel, has uniform control flow (no divergence), and coalesces the B read (consecutive col → consecutive addresses). And yet it will report a GFLOP/s number that’s a small single-digit percentage of the H100’s ~67 TFLOP/s FP32 peak.⁵ Why?

The access pattern — why it’s slow

B (K×N)

col j

C (M×N)

C[i][j]

One thread for C[i][j] reads ALL of row i of A (K elems) and ALL of col j of B (K elems) from HBM.
Every one of the N threads in row i re-reads the SAME row i of A. Every one of the M threads in col j
re-reads the SAME col j of B. Each A and B element is fetched from HBM ~N and ~M times → massive redundancy.

Figure 2 — Naive matmul’s global-memory access pattern. The thread owning C[i][j] streams a full row of A and a full column of B straight from HBM. Nothing is shared between threads: the N threads computing row i of C each independently re-read the identical row i of A, and the M threads computing column j each re-read the identical column j of B. Every input element is pulled from HBM roughly N (or M) times instead of once. Lesson 05 fixes exactly this by staging tiles in shared memory so a loaded element is reused by a whole block.

Compute the arithmetic intensity — and place it on the roofline

Let’s count, following Lesson 03’s method. Take $M = N = K = n$ for simplicity, FP32 (4 bytes).

FLOPs: each output element does $n$ multiply-adds $= 2 n$ FLOPs; there are $n^{2}$ outputs:
$FLOPs = 2 n^{3} .$

HBM bytes moved (naive, no reuse): each of the $n^{2}$ threads reads a full row of A ( $n$ elements) and a full column of B ( $n$ elements) $= 2 n$ elements $= 8 n$ bytes, plus one 4-byte write. Across all threads:
$bytes = n^{2} \cdot (8 n + 4) \approx 8 n^{3} .$

Arithmetic intensity:

I_{naive} = \frac{2 n ^{3}}{8 n ^{3}} = \frac{1}{4} = 0.25 \frac{FLOP}{byte} .

It’s constant — independent of $n$ . Now the roofline. Naive matmul runs on the FP32 CUDA cores (no Tensor Cores here), so the relevant peak is ~67 TFLOP/s and HBM is ~3.35 TB/s, giving an FP32 ridge point:⁵

I_{FP32}^{*} = \frac{67 \times 1 0 ^{12}}{3.35 \times 1 0 ^{12}} \approx 20 \frac{FLOP}{byte} .

Since $I_{naive} = 0.25 ≪ 20 = I^{*}$ , the kernel sits far out in the memory-bound region — by nearly two orders of magnitude. Attainable performance is capped at:

P \leq I \times β = 0.25 \times 3.35 TB/s \approx 0.84 TFLOP/s,

i.e. ~1.3% of the FP32 peak, before you even account for uncoalesced A reads and latency. The FP32 ALUs (and the Tensor Cores entirely) sit idle waiting on HBM. The kernel is not compute-bound; it’s starving.

The pathology is redundant HBM traffic from zero data reuse. Contrast with the potential intensity of matmul if you read each input only once: $bytes \approx 3 n^{2} \cdot 4 = 12 n^{2}$ , giving $I = 2 n^{3} /12 n^{2} = n /6$ — which grows with $n$ and, for $n = 1024$ , is ~170 FLOP/byte, firmly compute-bound. The naive kernel throws away that entire $O (n)$ reuse factor by re-fetching every element $O (n)$ times.

Quiz: Why is the naive matmul memory-bound even though matmul is the canonical compute-bound workload (Lesson 03)?

Answer

Matmul has potentially high arithmetic intensity — $O (n^{3})$ FLOPs over only $O (n^{2})$ distinct data — but that only materializes if each input element, once fetched from HBM, is reused across many output computations. The naive one-thread-per-output kernel does no reuse: the thread for $C_{ij}$ fetches row $i$ of A and column $j$ of B directly from HBM, and the $n$ different threads along a row of C each independently re-fetch the same row of A (and similarly for columns of B). So each element is read from HBM $\sim n$ times instead of once, inflating traffic from $O (n^{2})$ to $O (n^{3})$ and collapsing intensity to a constant $I = 0.25$ FLOP/byte — well below the FP32 ridge (~20). The FLOPs are inherent to the algorithm; the bytes are an artifact of the implementation. Fixing it means making threads cooperate to reuse loaded data — which requires shared memory and tiling (Lesson 05). (In practice the L2 cache recovers some of this reuse, so the measured kernel does better than the 0.25 back-of-envelope suggests — but you can’t rely on L2 at scale, and the kernel still leaves ~90%+ of peak on the table.)

This is the entire premise of Lesson 05

The naive kernel’s flaw is a memory-hierarchy flaw, not a math flaw. The fix — tiling through shared memory — loads a tile of A and a tile of B from HBM into SMEM once, then has every thread in the block reuse those tiles for a whole tile of outputs. That raises $I$ from a constant $0.25$ to $O (tile width)$ , marching the kernel right across the ridge into the compute-bound region where the ALUs (and eventually Tensor Cores) can actually be saturated. Everything you learned about SMEM, coalescing, and bank conflicts in Lesson 02 exists to make that tiled kernel fast.

10. Compiling with `nvcc`

nvcc is NVIDIA’s compiler driver. It splits your .cu file: host code goes to your system compiler (gcc/clang/MSVC), device code is compiled to PTX (a virtual ISA) and then to SASS (the actual machine code for a specific architecture).⁶

# basic
nvcc vadd.cu -o vadd
 
# optimized, targeting a specific architecture
nvcc -O3 -arch=sm_90 vadd.cu -o vadd     # sm_90 = Hopper (H100/H200)
#   sm_80 = A100 (Ampere) · sm_86 = RTX 30xx · sm_89 = L40/RTX 40xx · sm_90a = Hopper w/ wgmma
 
# see per-kernel register & shared-memory usage (occupancy input, Lesson 01 §7)
nvcc -O3 -arch=sm_90 --resource-usage vadd.cu -o vadd
#   or: nvcc ... -Xptxas -v   -> "Used N registers, M bytes smem" per kernel
 
# build a "fat binary" that runs on multiple GPUs (embeds SASS for several + PTX fallback)
nvcc -O3 -gencode arch=compute_80,code=sm_80 \
         -gencode arch=compute_90,code=sm_90 \
         -gencode arch=compute_90,code=compute_90 vadd.cu -o vadd

Key flags to know:

-arch=sm_XX — the compute capability of your target GPU. If you omit it you get an old default that won’t use newer features (and may not run at all on new hardware without JIT). Always set it.
--resource-usage / -Xptxas -v — prints registers/thread and shared memory/block per kernel. This is the input to the occupancy arithmetic from Lesson 01 §7 — check it whenever a kernel is slower than expected (look for register spills too).
-lineinfo — embed source line mapping so Nsight Compute and compute-sanitizer can point at your source lines. Cheaper than full -G.
-G — full device debug info (for cuda-gdb). Disables optimizations — never benchmark a -G build.

PTX vs SASS, and JIT

nvcc can embed SASS (final machine code, tied to one sm_XX) and/or PTX (forward-compatible virtual ISA). If a binary runs on a GPU it has no matching SASS for, the driver JIT-compiles the embedded PTX at load time (the first-launch cost you skip with a warm-up in timing).⁶ Shipping PTX (code=compute_90) buys forward compatibility to future GPUs; shipping SASS (code=sm_90) skips JIT for known ones. Fat binaries include both for a set of targets.

11. Practice & resources

Hands-on tasks

Write and run vector add yourself. Type out vadd.cu from §6 (don’t copy-paste — the muscle memory is the point), compile with nvcc -O3 -arch=sm_XX, and confirm max error = 0. Then add CUDA-event timing (§8) and print achieved GB/s. Compare to your GPU’s spec HBM bandwidth — you should hit a large fraction of it. Sweep threadsPerBlock ∈ {32, 64, 128, 256, 512, 1024} and note the (small) effect.
Convert it to a grid-stride loop. Rewrite the kernel as in §7, launch blocks = 32 * numSMs regardless of n, and verify correctness for n that is not a multiple of the block size (e.g. n = 1000003). Then launch it with <<<1, 1>>> and confirm it still produces the right answer (just slowly) — proof the loop decouples correctness from grid size.
Write and run naive matmul. Type out matmul_naive.cu from §9. Verify C[0] == K. Time it, compute achieved GFLOP/s, and divide by ~67 TFLOP/s (H100 FP32 peak) to get your % of peak — you should see a low single-digit-to-teens percentage. This number is your baseline for Lesson 05; write it down.
Profile the matmul. Run ncu --set full ./matmul_naive (Nsight Compute). Read the Memory Workload Analysis: check DRAM throughput %, global-load efficiency, and L2 hit rate. Look at the built-in roofline chart — the kernel’s dot should sit on the memory roof, confirming the §9 analysis empirically. Also run compute-sanitizer ./matmul_naive to confirm it’s memory-clean.
Where to run if you have no local GPU:
- Google Colab — free T4/occasionally better; Runtime → Change runtime type → GPU. Use !nvcc file.cu -o file && ./file in a cell, or the %%writefile magic to author .cu files. Zero setup.
- LeetGPU — browser-based CUDA playground with a growing set of kernel challenges (vector add, matmul, reductions). Great for drilling the exact kernels in this lesson with instant feedback.
- Lightning AI Studios / Modal — cheap on-demand A100/H100 for when you want real datacenter-GPU numbers (the 67 TFLOP/s and 3.35 TB/s figures only mean something on the actual hardware). Modal’s @app.function(gpu="H100") + a subprocess nvcc call is a clean way to run a .cu file on an H100 for cents.

Practice problem 1: Achieved bandwidth vs. arithmetic-intensity ceiling

Problem. You run vadd on an H100 (HBM ≈ 3.35 TB/s) with n = 2^26 floats and CUDA events report 0.28 ms for the kernel. (a) What effective bandwidth did you achieve, and what fraction of peak is that? (b) From the roofline, what’s the maximum FLOP/s this kernel could ever reach, and why is it so far below the 67 TFLOP/s FP32 peak?

Worked solution.

(a) Vector add touches 3 arrays (read A, read B, write C): bytes $= 3 \cdot 2^{26} \cdot 4 = 3 \cdot 268, 435, 456 \approx 8.05 \times 1 0^{8}$ B. Effective BW $= 8.05 \times 1 0^{8} /0.28 e- 3 \approx 2.87 \times 1 0^{12}$ B/s $\approx$ 2.87 TB/s, i.e. $2.87/3.35 \approx$ 86% of peak HBM. That’s a healthy memory-bound kernel — it’s doing about as well as a bandwidth-bound op can.

(b) Arithmetic intensity $I = \frac{1 FLOP}{12 bytes} \approx 0.083$ FLOP/byte. Roofline ceiling $= I \times β = 0.083 \times 3.35 e 12 \approx 2.8 \times 1 0^{11}$ FLOP/s $=$ ~0.28 TFLOP/s — about 0.4% of the 67 TFLOP/s FP32 peak. The kernel is doing almost no arithmetic per byte, so no amount of compute helps: it is fundamentally bandwidth-bound, and 86% of HBM is the win. This is the roofline telling you “stop optimizing compute; you’re already near the only ceiling that matters here.”

Practice problem 2: How much does tiling need to raise intensity?

Problem. For the naive FP32 matmul ( $M = N = K = 1024$ ) with $I_{naive} = 0.25$ FLOP/byte on an H100 (FP32 ridge $I^{*} \approx 20$ ), (a) by what factor must a tiled kernel raise arithmetic intensity to reach the compute-bound region? (b) A tiled kernel with tile width $T$ loads each element from HBM once per tile-row/col instead of once per output, giving $I \approx T /4 \cdot (something)$ — more precisely, staging $T \times T$ tiles makes each HBM-loaded element serve $T$ MACs, so $I \approx T$ ·(const). Roughly what tile width crosses the ridge, and what on-chip resource caps $T$ ?

Worked solution.

(a) From $0.25$ to $> 20$ is a ≥80× increase in arithmetic intensity. That’s the size of the gap tiling must close — a good gut-check on why the naive kernel is so far from peak (~1–2%).

(b) In classic tiled GEMM (Lesson 05), loading a $T \times T$ tile of A and of B into shared memory lets each loaded element participate in $T$ multiply-adds before being evicted, so HBM traffic drops by ~ $T \times$ and $I$ rises to roughly $I \approx \frac{2 n ^{3}}{2 n ^{3} / T \cdot 4} = T /4 \cdot 2 = T /2$ FLOP/byte (order-of-magnitude; the exact constant depends on layout). To exceed $I^{*} \approx 20$ you need $T ≳ 40$ -ish, but even a modest $T = 32$ gets you to ~16 FLOP/byte — right at the ridge and a ~64× improvement over naive, which in practice takes the kernel from a few percent to well over half of peak. What caps $T$ : shared-memory capacity per SM (two $T \times T$ FP32 tiles = $2 T^{2} \cdot 4$ bytes must fit in the ≤228 KB SMEM budget while leaving room for enough resident blocks)⁵ and the register file (each thread accumulates a sub-tile of outputs in registers). This is the exact occupancy-vs-working-set trade-off from Lesson 01 §7 — and the whole subject of the next lesson.

Real resources

📖 PMPP — Programming Massively Parallel Processors, 4th ed. (Hwu, Kirk, El Hajj). Ch. 2 (heterogeneous data-parallel computing — vector add, exactly this lesson’s skeleton) and Ch. 3 (multidimensional grids — the 2D matmul indexing). The canonical text; read these two chapters alongside this lesson.
📘 NVIDIA CUDA C++ Programming Guide — authoritative reference. Start with Programming Model (kernels, thread hierarchy, memory hierarchy) and Programming Interface (cudaMalloc/cudaMemcpy, launch, error handling, Unified Memory). https://docs.nvidia.com/cuda/cuda-c-programming-guide/
✍️ “An Even Easier Introduction to CUDA” (Mark Harris, NVIDIA Developer Blog) — the gentlest possible on-ramp; builds vector add with cudaMallocManaged and the grid-stride loop step by step, with timing. The perfect companion to §6–§8. https://developer.nvidia.com/blog/even-easier-introduction-cuda/
✍️ “How to Implement Performance Metrics in CUDA C/C++” and “How to Optimize Data Transfers” (NVIDIA blog) — the correct way to time kernels with CUDA events and reason about H2D/D2H cost. Reinforces §8.
🎥 GPU MODE (formerly CUDA MODE) lectures — the early lectures build exactly these first kernels hands-on and then march straight into optimizing matmul; the ideal next step after this lesson. https://github.com/gpu-mode/lectures
🔧 compute-sanitizer & Nsight Compute (ncu) docs — the two tools that turn “it’s wrong / it’s slow” into a specific diagnosis. Make both a reflex.

12. What’s next

Three threads, in order:

05-optimizing-gemm (next lesson) — take the naive matmul from §9 and make it fast: shared-memory tiling, register blocking, coalesced tile loads, and bank-conflict-free layouts (Lesson 02 §5–§6), marching it across the roofline ridge from ~1% toward >80% of peak. This lesson is the “before”; that one is the “after.”
06-memory-wall-and-flashattention (upcoming) — the same tiling/fusion idea applied to attention, the flagship IO-aware kernel.
Back to the diagnostic lens: 03-performance-modeling-roofline — re-read §9’s intensity calculation with the roofline math fresh; the naive matmul is the cleanest possible worked example of “compute the FLOPs and bytes, take the ratio, place it on the roofline, know your ceiling before you optimize.”

Topic hub: gpu-systems-for-llms | Lesson 04 of the GPU-systems curriculum | Filed: 2026-09-02
Related concepts: cuda-execution-model · gemm · memory-coalescing

References

NVIDIA CUDA C++ Programming Guide — authoritative reference for the execution-space qualifiers (__global__/__device__/__host__), the thread hierarchy and per-block thread limit (max 1024 threads/block; max grid dimension $2^{31} - 1$ blocks in x), the memory-management API (cudaMalloc/cudaMemcpy/cudaFree), Unified Memory / on-demand page migration (cudaMallocManaged), asynchronous kernel launches and error handling, and CUDA events for device-side timing. https://docs.nvidia.com/cuda/cuda-c-programming-guide/ ↩ ↩² ↩³ ↩⁴ ↩⁵ ↩⁶
Wen-mei W. Hwu, David B. Kirk & Izzat El Hajj, Programming Massively Parallel Processors: A Hands-on Approach, 4th ed. (Morgan Kaufmann, 2022). Ch. 2 (heterogeneous data-parallel computing — the allocate/copy/launch/copy-back vector-add skeleton) and Ch. 3 (multidimensional grids — the 2D one-thread-per-output-element matmul indexing). ↩ ↩²
Mark Harris, “CUDA Pro Tip: Write Flexible Kernels with Grid-Stride Loops,” NVIDIA Developer Blog (2013). Introduces the grid-stride loop idiom for decoupling grid size from problem size while preserving coalescing. https://developer.nvidia.com/blog/cuda-pro-tip-write-flexible-kernels-grid-stride-loops/ ↩
NVIDIA Compute Sanitizer documentation — the functional-correctness tool (formerly cuda-memcheck) that detects out-of-bounds/misaligned accesses, races, and uninitialized reads in kernels. https://docs.nvidia.com/compute-sanitizer/ ↩
NVIDIA H100 Tensor Core GPU Architecture Whitepaper (v1.01) and the H100 datasheet. H100 SXM5: FP32 (CUDA cores) ~67 TFLOP/s peak; 80 GB HBM3 at 3.35 TB/s; up to 228 KB configurable shared memory per SM. https://resources.nvidia.com/en-us-tensor-core and https://developer.nvidia.com/blog/nvidia-hopper-architecture-in-depth/ ↩ ↩² ↩³ ↩⁴
NVIDIA CUDA Compiler Driver NVCC documentation — nvcc splits host/device code, compiles device code to PTX (virtual ISA) and SASS (architecture-specific machine code), embeds fat binaries, and relies on driver JIT of PTX for forward compatibility. https://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/ ↩ ↩²

Gunjan Dhanuka — Learning Notes

Explorer

Learn: Your First CUDA Kernels

Learn: Your First CUDA Kernels

1. Learning map

2. Why this matters

3. The CUDA programming model: host vs device

Function qualifiers

4. Launch configuration and the global thread index

The thread-index → data-element mapping (vector add)

The boundary guard

5. Memory management

Unified Memory: `cudaMallocManaged`

6. Vector add — end to end

7. Grid-stride loops

8. Error checking, async launches, and synchronization

The `CUDA_CHECK` macro

Checking a kernel launch — the two-step

Timing with CUDA events

9. Naive matmul — correct, and instructively slow

The access pattern — why it’s slow

Compute the arithmetic intensity — and place it on the roofline

10. Compiling with `nvcc`

11. Practice & resources

Hands-on tasks

Real resources

12. What’s next

References

Graph View

Table of Contents

Backlinks

Gunjan Dhanuka — Learning Notes

Explorer

Learn: Your First CUDA Kernels

Learn: Your First CUDA Kernels

1. Learning map

2. Why this matters

3. The CUDA programming model: host vs device

Function qualifiers

4. Launch configuration and the global thread index

The thread-index → data-element mapping (vector add)

The boundary guard

5. Memory management

Unified Memory: cudaMallocManaged

6. Vector add — end to end

7. Grid-stride loops

8. Error checking, async launches, and synchronization

The CUDA_CHECK macro

Checking a kernel launch — the two-step

Timing with CUDA events

9. Naive matmul — correct, and instructively slow

The access pattern — why it’s slow

Compute the arithmetic intensity — and place it on the roofline

10. Compiling with nvcc

11. Practice & resources

Hands-on tasks

Real resources

12. What’s next

References

Footnotes

Graph View

Table of Contents

Backlinks

Unified Memory: `cudaMallocManaged`

The `CUDA_CHECK` macro

10. Compiling with `nvcc`