100 DAYS OF GPU Course preview
Day 43Module 5
in-technical-review

Optimizing CUDA matmul, steps 1 to 3

Day 16 put a 16 by 16 tile under a matrix multiply, cut global loads per thread by exactly sixteen times, and recorded this result:

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

Sixteen times fewer loads gave 1.53 times the speed. Two facts explain the difference, and a profiler measures both.

The first is that day 16's naive kernel was not the naive kernel. It was already coalesced, which is step 2 of this three-step sequence. The 1.53x result measures only the final change.

The second fact is that at 256 by 256 the whole problem fit in L2, so the repeated reads tiling removes were never crossing to memory anyway.

This page builds all three steps, profiles each one, and grades them on four ratios instead of on milliseconds.

What decides a warp's 32 addresses in a matmul

One thread owns one element of C. It walks a row of A and a column of B, multiplies, adds, and stores once. The arithmetic is fixed.

The programmer chooses which built-in index becomes the row and which becomes the column, and that choice decides everything a warp asks global memory for.

A 16 by 16 block linearises as threadIdx.x + 16 * threadIdx.y, so one warp is sixteen consecutive x and two consecutive y. Global memory is served in 32-byte sectors, and a sector holds eight floats. Now count the sectors one warp's load request covers, for an N by N matrix of floats.

Step 1, threadIdx.x chooses the row. The warp sits on sixteen consecutive rows and two consecutive columns. a[row * n + p] is sixteen addresses 4N bytes apart, so sixteen sectors, and fifteen sixteenths of every one of them is bytes nobody asked for. b[p * n + col] is two adjacent floats inside a single sector.

Seventeen sectors for the pair of requests, so 8.5 per request.

Step 2, threadIdx.x chooses the column. The warp sits on two rows and sixteen consecutive columns. a[row * n + p] is now two addresses, one sector each. b[p * n + col] is sixteen consecutive floats starting on a 64-byte boundary, which is two sectors, and both halves of the warp read the same two.

Four sectors for the pair, so 2.0 per request.

That is a factor of 4.25 in sectors moved, from two lines. The request count does not move at all, which is the same result day 11 measured on a copy kernel: the access pattern cannot change how many instructions issue, only what each one costs.

Step 3 stages both inputs through shared memory. Each thread now issues two loads per tile step instead of two per k step, so 2 * ceil(N / 16) instead of 2N. Both tile loads are two rows of sixteen consecutive floats, which is four sectors each, so 4.0 per request. Higher than step 2, on one sixteenth as many requests.

The step that usually gets skipped

Many introductions compare an initial kernel with a shared-memory kernel. They often omit whether the first kernel already coalesces its memory access.

const size_t row = blockIdx.x * blockDim.x + threadIdx.x;  // step 1
const size_t col = blockIdx.x * blockDim.x + threadIdx.x;  // step 2

Both versions give the correct answer. One moves 4.25 times as many sectors as the other. A source review may miss the difference because only the index assignment changes.

This is why CUDA-CODE-STYLE.md pins row to y and col to x for the whole course. The convention prevents this error.

You still need the comparison to tell whether tiling helped through reuse or whether it also fixed the baseline's access pattern.

Day 16 used step 2 as its baseline. Its 1.53x result measures the change from step 2 to step 3.

One variable per step

Full program in matmul_steps.cu, under code/day43-matmul-1/. Four choices keep the comparison controlled.

One thing moves per step. Step 1 and step 2 differ by two lines. The launch configuration, the block shape, the arithmetic, the guard and the answer are identical, and because both cases are square the grid is the same either way. Step 3 adds the shared tile and changes nothing else.

Steps 2 and 3 are day 16's kernels, with the three dimensions collapsed to one because every case here is square. The N = 512 rows of this program and day 16's 512 rows are two independent measurements of one pair. If they disagree, one of the two programs is wrong, and that is worth more than either number.

Two sizes, checked against the GPU's L2 size. Three 512 by 512 float matrices use 3 MiB; three 1024 by 1024 matrices use 12 MiB. A static_assert checks these byte counts.

The program reads the L2 size from the device and notes whether each case fits. On a GPU with at least 12 MiB of L2, the second case is no longer an out-of-cache test.

Events, a warm-up per kernel, and all three checked before any is timed. That last part is not tidiness. It means the first three kernel launches of the process are one of each step, which is what lets ncu --launch-count 3 collect exactly the set this page quotes.

Step 1 is the whole of the difference:

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

Step 2 is that kernel with these two lines in place of its first two:

    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;

Note. The sectors-per-request column the program prints is arithmetic from the block shape, not a measurement. It is printed so the Nsight Compute report has something to be checked against, the same way day 16 printed loads per thread. If the report disagrees with it, the model in the section above is wrong and the page changes, not the report.

Results

Measured. Tesla T4, driver 595.84, CUDA 12.6 (V12.6.85), built with nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo, profiled with Nsight Compute 2024.3.2 under sudo at --clock-control base. Captured 2026-09-01 on the project's verification node; the transcript is in the page's evidence file and both reports ship at content/profiles/optimizing-matmul-1/.

Re-verified on 2026-09-02 with CUDA 13.0 (V13.0.88), driver 580.173.02, on the same Tesla T4. All three kernels passed at both sizes and the CUDA event ratios stayed in the same bands: 8.38x end to end at N = 512 and 7.87x at N = 1024. The fresh Nsight Systems capture is listed in evidence; the counter table below remains the CUDA 12.6 Nsight Compute capture.

The timed table, CUDA events, mean of 10 runs after 3 warm-ups:

   N  kernel                 time (ms)     GFLOP/s   loads/thr  sect/req
----  -------------------  -----------  ----------  ----------  --------
 512  matmulStridedRows          5.381        49.9        1024       8.5
 512  matmulCoalescedCols        1.024       262.2        1024       2.0
 512  matmulTiled                0.645       416.5          64       4.0
      speedup  step1->step2  5.26x  step2->step3  1.59x  step1->step3  8.35x
1024  matmulStridedRows         22.952        93.6        2048       8.5
1024  matmulCoalescedCols        4.003       536.5        2048       2.0
1024  matmulTiled                2.735       785.2         128       4.0
      speedup  step1->step2  5.73x  step2->step3  1.46x  step1->step3  8.39x

And the report, N = 512 launches (the N = 1024 rows sit in the second shipped report):

Metric step 1 step 2 step 3
sm__throughput.avg.pct_of_peak_sustained_elapsed 11.58 60.97 72.63
gpu__compute_memory_throughput.avg.pct_of_peak_sustained_elapsed 49.19 60.97 72.63
sm__warps_active.avg.pct_of_peak_sustained_active 94.43 95.33 94.87
l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum 8,388,608 8,388,608 524,288
l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum 71,303,163 16,777,072 2,084,470
lts__t_sector_hit_rate.pct 96.96 97.70 97.30
dram__bytes_read.sum + dram__bytes_write.sum 7.82 MB 4.71 MB 4.60 MB

The prediction, against the card

  1. Matched. Sectors per request: 8.50, 2.00, 3.98 at N = 512 and 8.50, 2.00, 3.97 at N = 1024. The tiled kernel's 3.98 is the sector counter reading a handful under its derived 4.0, the same sliver day 42 saw, not a geometry error.
  2. Matched exactly. 8,388,608 requests for step 1 and step 2 at N = 512, 524,288 for step 3. The feared 128-bit unrolling did not happen: nvcc could not prove alignment for a runtime n, exactly as the caveat guessed.
  3. Matched. Step 1 to step 2 is 5.26x at 512 and 5.73x at 1024; step 2 to step 3 is 1.59x and 1.46x. Neither matches its sector ratio (4.25 and 8): time is not bytes.
  4. Did not match. Tiling won by less at N = 1024, 1.46x against 1.59x, not more. The prediction's mechanism was real but aimed at the wrong step: at 1024 the L2 hit rate for step 2 fell from 97.70 to 83.42 percent, so the coalesced kernel got slower per element and step 2's own speedup over step 1 grew instead (5.26x to 5.73x). What out-of-L2 bought was taken by step 2, not left for step 3.
  5. Matched. Achieved occupancy spans 94.43 to 95.33 percent at N = 512 while the times span 8.35x. Step 1 keeps its warps resident and stalled, latency hiding failing while the occupancy number looks healthy; the healthiest-looking number in the table belongs to no ranking of speed.

Do not compare your milliseconds with these results. Absolute times vary by GPU, so compare these four ratios: speed-of-light percentage, achieved occupancy, sectors per request, and L2 hit rate.

Run it yourself

Run the timing test on a CUDA GPU. Add -lineinfo because this lesson also uses a profiler, and change -arch for your GPU if needed:

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

Compiler Explorer is the wrong tool for this page. The harness runs two CPU references, checks three kernels at two sizes, and issues more than sixty timed launches; that is a harness plus data, not the single-buffer, single-kernel shape one editor pane and a 20 second cap were built for.

Hardware. Nsight Compute needs counter permission. Plain ncu was blocked on this project's node because RmProfilingAdminOnly was 1. The current free-Colab case remains untested. Do the whole exercise from the reports shipped with this page, at content/profiles/optimizing-matmul-1/: optimizing-matmul-1-512.ncu-rep, optimizing-matmul-1-1024.ncu-rep, their .details.txt and .details.csv exports, and optimizing-matmul-1.nsys-rep. Nsight Systems needs no counters and runs anywhere. With no GPU at all, start at /setup/learn-cuda-without-a-gpu.

Exercise

Starting from matmulStridedRows, write step 2 and then step 3, and report the two speedups your card gives you alongside the sectors per request each kernel reaches in the shipped report.

Time: 30 to 45 minutes. Submit: your matmul_steps.cu, the two speedups, the three sectors-per-request figures, and one sentence saying which step you would have skipped and what it would have cost you.

Check: the harness runs both of your kernels against the CPU reference at N = 512 and N = 1024 and, on a mismatch, prints the kernel, the case, and the first row and column that disagree. It then measures step 1 as its own baseline on your card, in the same process, seconds before it times your two, and requires step 2 to beat step 1 at both sizes and step 3 to beat step 2 at N = 1024, the size that does not fit L2.

At N = 512 the step 3 ordering is reported, not gated: 3 MiB of matrices inside a 4 MiB L2 is exactly the regime where day 34 measured a tile losing to the cache it duplicates, and this page's own prediction 4 allows the same outcome here. It reports every ratio either way, because the ratios are the answer and passing is only the gate.

It also runs the memory-traffic audit and holds step 3 to 2 * ceil(N / 16) + 8 global reads per thread, which is 72 at N = 512, so a correct step 2 submitted as step 3 passes correctness and still fails the day, and the report says so.

The tile stays at 16 by 16. A static_assert pins it, because the sectors-per-request figures the program prints are derived for that block and changing it means re-deriving all three. Day 44 is where the tile moves.

Hint 1

The two kernels launch identically, compute identically and agree at every element. Look only at which built-in index ends up in the first subscript of a, and ask what the 32 lanes of one warp are doing at a single instant.

Hint 2

One warp is threadIdx.x + 16 * threadIdx.y for x in 0 to 15 and y in 0 to 1. A 32-byte sector holds eight floats.

Count the distinct sectors A's load touches under each of the two assignments. One answer is sixteen and the other is two.

Solution

Step 2 is step 1 with the two index lines swapped, so col comes from x and row comes from y. Nothing else changes, and the snippet above is the whole diff. Step 3 is day 16's matmulTiled: two __shared__ arrays, a loop over ceil(N / 16) tile steps with a __syncthreads() on each side of the inner product, and a load that writes 0.0f rather than skipping when an element is out of range.

The three ratios to report are step 1 to step 2, step 2 to step 3, and the same pair again at the larger size. Expect the first to be the larger of the two, and expect the second to grow when the working set stops fitting your card's L2. Your absolute times will not match anyone else's; the ordering and the ratios will.

Coalescing and reuse reduce different parts of memory traffic. Step 2 fixed the shape of each request and left the count alone. Step 3 divided the count by the tile width and made each surviving request larger.

A kernel may need both changes. Compare sectors per request with request count to identify which change is missing.

Pitfalls

Your tiled kernel is only about 1.5x faster and you expected 16x. Two causes can combine. The kernel you called naive was probably already coalesced, so you measured the last step on its own; and if your matrices fit your card's L2, the repeated reads tiling removes were cache hits rather than memory traffic.

Measure step 1 as well, and run a size past your L2. Day 30 is where the ceiling comes from and day 42 is where the memory chart shows it.

ncu prints nothing but a permission error. The string is ERR_NVGPUCTRPERM, in full: The user running <tool_name/application_name> does not have permission to access NVIDIA GPU Performance Counters or the Hardware Event System on the target device. The cause is the stock driver default RmProfilingAdminOnly: 1 in /proc/driver/nvidia/params, so sudo ncu works where you have root. The current free-Colab case remains untested; use the shipped reports when local counter permission is unavailable.

You compare your milliseconds against the ones on this page. A T4 and an H100 differ by more than an order of magnitude on every absolute counter here, so a target time would fail an H100 owner for owning an H100. Compare the four ratios instead.

Achieved occupancy went up and the kernel got slower. Occupancy counts warps resident on the SM, not warps issuing instructions, and a kernel stalled on memory keeps every warp it has. It is a means to latency hiding, never a score. Day 45 makes that trade the whole subject.

Sectors per request went up when you tiled, so you undid the tiling. It should go up. Total sectors is requests times sectors per request, and tiling divides the first by sixteen while doubling the second. Read the two numbers together or neither means anything.

Day 11 is where that pair is introduced.

You profiled a warm-up launch instead of the kernel you meant. --launch-count counts from the first launch of the process, so a program that times each kernel as it reaches it hands you thirteen launches of the first one. This program checks all three before timing any, which is why --launch-count 3 is enough; anything else needs --launch-skip.

Go deeper

Next

Day 44 applies steps 4 to 6. Each thread computes several outputs from registers, and the loads become wider. The target is a fraction of cuBLAS rather than an absolute time.

Day 45 tunes one kernel to the same speed at two occupancy values. Day 49 places all six kernels on a roofline built from measured compute and bandwidth limits, then uses arithmetic intensity to classify each memory-bound result.