100 DAYS OF GPU Course preview
Day 16Module 2
in-technical-review

CUDA tiled matrix multiplication, measured

A naive matrix multiply at 512 by 512 makes each thread read 1,024 floats out of memory to produce one output element. The tiled version reads 64. Same arithmetic, same answer, one sixteenth of the loads, and the sixteen is the tile width.

That ratio explains why tiling exists. It is also why so much tiled matmul code only works when the tile width divides the matrix. The guard for the leftover edge is a few characters, and it is not the guard people write first: write the obvious one and your kernel passes at 512, passes at 1024, and multiplies stale shared memory into the answer at 611.

This page builds both kernels, counts what each one loads, and runs them at three sizes, one of which divides by 16 on no axis.

Where a naive matmul spends its loads

C = A * B with A of shape M by K and B of shape K by N. One output element is a dot product: walk a row of A and a column of B, multiply, add. The obvious CUDA version gives one thread one output element, so that thread does 2K reads and 2K floating point operations, then writes once.

Now count the whole grid. Row 0 of A is read by every thread that owns an element of row 0 of C, and there are N of them. Column 0 of B is read by M threads.

Every input element crosses from global memory into the chip O(N) times, and none of those threads knows the others exist.

Put a number on it the way day 9 did. Each thread issues 2K loads, 4 bytes each, and does 2K operations, so it does 0.25 floating point operations for every byte it asks for. That figure is arithmetic intensity, and 0.25 is low enough that no card in this course's hardware matrix can hide the traffic behind the arithmetic.

The kernel is memory-bound and stays memory-bound however you tune the launch.

The threads that read the same value are not scattered across the machine. The 256 threads of one thread block covering a 16 by 16 patch of C between them need exactly 16 rows of A and 16 columns of B, and they need them at the same time.

They share a scratchpad, shared memory, which is on the SM rather than across the bus. Load a 16 by 16 square of A and a 16 by 16 square of B into it once, have every thread do its 16 multiply-adds out of the squares, then move on to the next square along K.

That is tiling. Each thread now issues 2 loads per tile step instead of 2 per k step, so 2 * ceil(K / 16) in total. The arithmetic did not change, so the intensity goes from 0.25 to 4 operations per byte: a factor of exactly the tile width.

It does not make the kernel compute-bound, it makes it 16 times less memory-bound, and a wider tile buys more of the same.

Why a tiled kernel cannot use the guard you know

Days 5 and 7 taught one shape for a bounds check: compute your index, test it, and put the work inside the test. That shape is right in every kernel so far, because none of them had a barrier.

A tiled kernel has two. Every thread has to finish writing its element of the tile before any thread reads the tile, and every thread has to finish reading the tile before any thread overwrites it on the next step. Both of those are __syncthreads().

A __syncthreads() that only part of a block reaches is undefined behaviour, not a slow path. The shape below is not an option:

if (row < m && col < n) {
    // loads, __syncthreads(), multiply-adds, __syncthreads()
}

At 512 by 512 nothing goes wrong, because every block is full and the test is true for all 256 threads. At 611 by 613 the blocks along the bottom and right edges are ragged, some of their threads fail the test, and those threads never arrive at the barrier their neighbours are waiting on.

The rule is the one day 14 states and this kernel obeys: guard the loads and the store, never the barrier. Every thread runs the tile loop and reaches both barriers. A thread whose element is off the edge of A still writes something into the tile, and what it writes is the second half of the problem.

Two kernels, three sizes, one of them ragged

Full program in code/day16-tiled-matmul/tiled_matmul.cu. The comparison follows three rules.

Neither kernel wins by computing less. Same 16 by 16 block, same grid rounded up on both axes, and both produce every element of C. Only the load count differs.

The two square cases do exactly 2 * M * N * K operations each; on the ragged case the tiled kernel does about 5 percent more, because its edge tiles are padded with zeros and it multiplies those too. The GFLOP/s column divides by 2 * M * N * K in every row, so that 5 percent is charged to the tiled kernel rather than credited to it.

Both are checked before either is timed, against a single-threaded CPU reference at every element and at every size. The inputs are small integers held as floats, so the largest dot product any case reaches is 6 * 617 = 3,702, well under the 2^24 where a float stops representing consecutive integers. Every intermediate on both processors is exact, and a mismatch can only be an indexing bug.

Three sizes, and the third is the point. 256 and 512 divide by 16 on every axis. 611 by 613 by 617 divides on none: 611 is 38 whole tiles and 3 rows, 613 is 38 and 5, 617 is 38 and 9.

The three remainders differ, so a guard that is right on one axis and wrong on another cannot pass by luck. A static_assert holds both facts, because they follow from constants alone.

Here is the naive kernel:

__global__ void matmulNaive(const float* __restrict__ a,
                            const float* __restrict__ b, float* __restrict__ c,
                            size_t m, size_t n, size_t k) {
    const size_t col =
        blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    const size_t row =
        blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
    if (row < m && col < n) {
        float acc = 0.0f;
        for (size_t p = 0; p < k; ++p) {
            acc += a[row * k + p] * b[p * n + col];
        }
        c[row * n + col] = acc;
    }
}

It guards the whole body, and it is allowed to, because it has no barrier. Its loads are not badly laid out either: consecutive threads take consecutive columns, so b[p * n + col] is coalesced. The cost here is repeats, not pattern.

And the tiled one:

__global__ void matmulTiled(const float* __restrict__ a,
                            const float* __restrict__ b, float* __restrict__ c,
                            size_t m, size_t n, size_t k) {
    __shared__ float tileA[kTileDim][kTileDim];
    __shared__ float tileB[kTileDim][kTileDim];

    const size_t col =
        blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    const size_t row =
        blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;

    // ceilDiv is a host function, so the same arithmetic is written out here.
    const size_t tiles = (k + kTileDim - 1) / kTileDim;

    float acc = 0.0f;
    for (size_t tileIdx = 0; tileIdx < tiles; ++tileIdx) {
        const size_t aCol = tileIdx * kTileDim + threadIdx.x;
        const size_t bRow = tileIdx * kTileDim + threadIdx.y;

        tileA[threadIdx.y][threadIdx.x] =
            (row < m && aCol < k) ? a[row * k + aCol] : 0.0f;
        tileB[threadIdx.y][threadIdx.x] =
            (bRow < k && col < n) ? b[bRow * n + col] : 0.0f;
        __syncthreads();

        for (int p = 0; p < kTileDim; ++p) {
            acc += tileA[threadIdx.y][p] * tileB[p][threadIdx.x];
        }
        __syncthreads();
    }

    if (row < m && col < n) {
        c[row * n + col] = acc;
    }
}

Those two ternaries took two attempts. The obvious way to guard a load is to guard it:

if (row < m && aCol < k) {
    tileA[threadIdx.y][threadIdx.x] = a[row * k + aCol];
}

No out-of-range read, every thread still reaches both barriers, and it is still wrong. The tile is reused on every step of the k loop, so skipping the write leaves the value the previous step put there, and the inner loop multiplies that into the accumulator.

On 611 by 613 by 617 the last k step covers columns 608 to 623 of a matrix with 617 of them, so seven of tileA's sixteen columns still hold the step before. Writing 0.0f is not defensive coding: zero is the identity for a sum, so the padded lanes contribute nothing.

This is the same error as day 11's benchmark bug, made while writing the lesson that warns about it.

Note. The loads/thread column the program prints is arithmetic, not a measurement: 2K and 2 * ceil(K / 16) from the case dimensions. It counts loads a thread issues, not bytes DRAM moved. The L1 and L2 caches already serve some of the naive kernel's repeats, so the measured speedup will be smaller than that ratio. How much smaller is the interesting number here.

Results

Originally measured on a Tesla T4 with driver 595.84 and CUDA 12.6. The exact shipped program was re-verified with driver 580.173.02 and CUDA 13.0 (V13.0.88) on 2026-09-02, with no correctness or meaningful behavioral drift. All three cases still passed, the load reductions remained 16.0x, 16.0x and 15.8x, and the speedups remained 1.53x, 1.59x and 1.29x.

The original 12.6 capture below and both CUDA 13 captures remain in evidence.

GPU: Tesla T4 (compute capability 7.5)
tile 16 x 16, 256 threads per block, 2048 bytes of shared memory per block

case             kernel    time (ms)     GFLOP/s  loads/thread
---------------- ------- ----------- ----------- -------------
256x256x256      naive         0.165       203.5           512
256x256x256      tiled         0.108       312.0            32
                 tiled is 1.53x faster and reads 16.0x fewer floats per thread
512x512x512      naive         0.982       273.4          1024
512x512x512      tiled         0.620       432.9            64
                 tiled is 1.58x faster and reads 16.0x fewer floats per thread
611x613x617      naive         1.759       262.8          1234
611x613x617      tiled         1.364       338.9            78
                 tiled is 1.29x faster and reads 15.8x fewer floats per thread

all 3 cases match the CPU reference at every element

The loads-per-thread column is the lesson; the time column is the consequence. A 16 by 16 tile means each thread reads 16 times fewer floats from global memory, and the measurement says 16.0x on the two square cases, exactly the tile width. That is not an approximation, it is what tiling does: each value loaded once into shared memory is used by all 16 threads in the row that need it.

The speedup is smaller than the load reduction, 1.53x rather than 16x, and the gap matters. Cutting global traffic by 16x only helps until something else binds, and at these sizes the kernel is no longer purely memory-bound.

The non-square case is the one worth keeping. 611 by 613 by 617 is divisible by no tile dimension, so every tile on the right and bottom edges is partial. It still matches the CPU reference at every element, because the load guards write zeros into the out-of-range part of the tile rather than skipping the load.

The tiled speedup drops to 1.29x, since the edge tiles do full work for partial output. A tiled matmul that only works on multiples of the tile width is the classic version of this kernel and it is wrong.

None of these numbers is close to cuBLAS, and they are not meant to be. Day 43 starts the next set of optimizations.

Run it yourself

A free Colab T4 or any card you own. The build line is the one in the repo's README:

nvcc -std=c++17 -O3 -arch=sm_75 -o tiled_matmul tiled_matmul.cu

There is no Compiler Explorer embed on this page. The slow part is not the GPU, it is the CPU reference: a single-threaded triple loop over 512 by 512 by 512 and again over 611 by 613 by 617, which sits too close to Compiler Explorer's 20 second run cap to pin a lesson to. Cut the case list to the smallest size and it fits an embed fine, as long as you target sm_75 or lower, because that runner is a Tesla T4 too.

Exercise

Write the tiled kernel yourself, from a starter that hands you the naive one, and make it right on a size the tile does not divide.

Time: 30 to 45 minutes. Submit: your tiled_matmul.cu, plus one sentence naming what your guard writes into the tile when the element is off the edge, and why.

Check: the harness runs the GEMM case ladder, which includes (m, n, k) = (611, 613, 617), and reports the smallest failing case rather than the first, so a kernel that is right on multiples is reported at 611 instead of 4096.

It also prints a memory-traffic score: global reads per thread, against a budget of 72 at tile 16 and N 512. The naive kernel reads 1,024, so a correct submission that never touched shared memory passes correctness and still fails the day, and the report says so.

The full contract is at /reference/harness.

Hint 1

Count how many times one element of A crosses into the chip while the naive kernel runs. Then ask what the other 255 threads of that block want at the same moment.

Hint 2

The block fills a tile, every thread reads sixteen values out of it, then the block fills the tile again. Two of those transitions need something to happen first. Which two, and what happens to a thread that ran ahead?

Solution

matmulTiled in tiled_matmul.cu is the reference version, and the diff against the naive kernel is three things: two __shared__ arrays, a loop over ceil(K / 16) tile steps with a __syncthreads() on each side of the inner product, and a load that writes 0.0f rather than skipping when the element is out of range. The guard covers the loads and the final store, never the barriers. On the T4 in the table above that version runs 611 by 613 by 617 in 1.364 ms against the naive kernel's 1.759, which is where your own ragged row should land relative to your own naive row.

Keep this rule: you did not cut the arithmetic. You used each fetched value for sixteen multiply-adds instead of one. Days 43 and 44 apply the same rule to registers and wider loads.

Pitfalls

Your kernel passes at 1024 and is wrong at 611. Almost always the tile load: guarded with an if and no else, so a lane that is off the edge leaves the previous step's value in shared memory. Write 0.0f in the else.

Pick a test size no tile divides so the test exposes the bug.

The kernel hangs, or its answer changes between runs, on the ragged size. A __syncthreads() inside a branch only part of the block takes. Guard the loads and the store, never the barrier.

compute-sanitizer --tool synccheck names it, and day 14 explains the rule.

You compare a tiled kernel that skips the edge against a naive one that does not. That benchmark is timing two different amounts of work, and the wrong one wins. Check the answer before you time it, at a size that is not a multiple.

You expect the speedup to match the ratio of loads. It will not, because the caches were already serving some of the naive kernel's repeats. Tiling turns reuse you were getting by luck into reuse you are getting by construction.

The reuse still works as the matrix grows past the cache. Day 42 puts the two memory charts side by side.

Your tiled kernel is correct and no faster. Check that consecutive threads still take consecutive columns. Index the tile as tileB[threadIdx.x][p] and the 32 lanes of a warp collide in shared memory instead of spreading across banks, which is a bank conflict and cancels the gain from tiling.

Day 15 measures it.

Go deeper

Next

Day 17 reads the ptxas report for this kernel, because registers are the other resource a tile competes for. Day 20 closes the module by tiling a 2D convolution, where the tile needs a halo and stops lining up with the outputs. Days 43 and 44 take this kernel through register tiling and wider loads.