Case study: why naive matmul is slow

01 · The benchmark

The kernel everyone writes first

Multiply two 4096×4096 FP32 matrices. The work is 2N³ FLOPs — two operations per multiply-accumulate, 4096³ multiply-accumulates ≈ 137 GFLOPs. At the 4070's ~29.1 TFLOPS peak, the floor on runtime is about 4.7 milliseconds. Hold that number; the naive kernel will not be near it.

// one thread per output element C[r][c]
__global__ void matmul_naive(const float* A, const float* B,
                                    float* C, int N){
    int r = blockIdx.y * blockDim.y + threadIdx.y;
    int c = blockIdx.x * blockDim.x + threadIdx.x;
    if (r < N && c < N){
        float acc = 0.0f;
        for (int k = 0; k < N; ++k)
            acc += A[r * N + k] * B[k * N + c];
        C[r * N + c] = acc;
    }
}

Sixteen million threads, each marching down a row of A and a column of B. It's correct, readable, and — as day 12 of a journey about slowness — educational.

02 · Walking it through days 7–11

First surprise: the memory access is fine

Run day 07's checklist on the inner loop. Within a warp, r is fixed and c varies by 1 per lane, so the B[k*N + c] loads across the warp are 32 consecutive floats — perfectly coalesced. The A[r*N + k] load is identical for all 32 lanes — a broadcast, which GPUs serve once. No uncoalesced disaster. The problem is upstream of the access pattern, and day 10's classification finds it:

Count the reuse. Every element of A is read by N different threads (all the threads in its row), and every element of B by N threads (all the threads in its column). Total DRAM traffic, if nothing caches it: 2N³ × 4 bytes — for N = 4096, that's 549 GB. Per run. The kernel moves its own weight in data more than 45 times over just to do one pass of arithmetic.

03 · Placing it on day 11's chart

0.25 FLOPs per byte has nowhere to go

The arithmetic intensity falls out of the same count: 2 FLOPs per multiply-accumulate, and 8 bytes of traffic per multiply-accumulate (two 4-byte loads, one 4-byte store). That's 0.25 FLOPs/byte — deeply memory-bound, a ceiling of 504 GB/s × 0.25 ≈ 126 GFLOPS, about 0.43% of peak. The runtime floor from bandwidth alone: 549 GB ÷ 504 GB/s ≈ 1.09 seconds — 230× slower than the compute floor we computed in section 01.

04 · Where the fix lives

The problem is a data-flow problem, not a code problem

Notice what isn't wrong. The threads are mapped fine (day 03/14 territory), the warp isn't diverging (day 05), occupancy will be healthy, and the loads coalesce (day 07). By every check from the first week, this kernel passes. It fails only on reuse — and reuse is fixed by restructuring, not by tweaking: stage tiles of A and B in shared memory (day 08) so each byte loaded is used by many threads before being discarded.

That restructure is day 15, and it will lift the ceiling 16× with one idea. Beyond that lie register blocking, double buffering, and tensor cores — the gap between our tiled version and cuBLAS's ~85–90% of peak. But the diagnosis you just watched — count traffic, compute intensity, place on the roofline, fix the binding constraint — is the whole method. Days 13 and 14 first need to make sure you can actually launch these things; then day 15 fixes this one properly.

Next: day 13 — the Writing CUDA phase begins. Your first kernel, end to end: the host side, the launch syntax, thread indexing, and the error checking everyone skips.