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

Optimizing CUDA matmul, steps 4 to 6

Day 30 measured the tiled matmul at 976.0 GFLOP/s, or 101.3 percent of its measured bandwidth limit. The same run measured an FMA peak of 8092.1 GFLOP/s.

The kernel reaches its bandwidth limit while using 12 percent of the measured FMA rate. This page changes the design in steps 4 to 6 and compares each step with cublasGemmEx using the same FP32 arithmetic.

Two ladders, one kernel

Step 3 gave every thread one element of C. Each step of its inner loop issues two shared memory loads and one FMA, and that ratio of 2.0 load instructions per multiply-add is the whole problem. Loads and FMAs compete for issue slots, so a kernel built at 2.0 spends most of its instruction stream fetching operands it uses exactly once.

The next step reuses values in registers, avoiding a shared-memory load for each use. Give a thread a column of 8 outputs (step 4) and one value of B fetched from the tile feeds 8 FMAs: 9 loads per 8 multiply-adds, 1.125.

Give it an 8 by 8 micro-tile (step 5) and 8 values of A plus 8 of B feed 64 FMAs: 0.25. Widen those 16 scalar loads into 4 float4 loads (step 6) and the instruction count falls to 0.0625 per FMA, with the values read unchanged.

Diagram: what one thread holds in registers. Three bands, one axis: outputs owned per thread, with the shared traffic each one needs per loop step drawn as arrows from the two tiles into the thread's registers. Band 1: step 3, one output. Two loads feed one FMA. "2.0 loads per FMA." Band 2: step 4, a column of eight. Nine loads feed eight FMAs. "1.125 loads per FMA." Band 3: step 5, an eight by eight micro-tile. Sixteen loads feed sixty-four FMAs. "0.25 loads per FMA." Alt text: One output per thread costs two shared loads per multiply-add. A column of eight costs about one. An eight by eight micro-tile costs a quarter: sixteen values held in registers feed sixty-four multiply-adds.

The block tile also grows. A 16 by 16 block tile moves 257 global floats per output at N = 2048; the 128 by 128 tile of steps 5 and 6 moves 33. That takes arithmetic intensity from 4.0 FLOP per byte to 31.0, and day 30 measured this card's ridge at 33.09.

The larger tile increases global-memory reuse, while the micro-tile increases register reuse. Both changes are needed to move the kernel away from a memory-bound limit.

Step 6 adds one mechanism. Its inner loop needs a thread's 8 values of A adjacent in shared memory, so the tile is stored transposed, and the fill that writes it scatters each float4's four lanes down a column. Two warp lanes always land 512 floats apart, and 512 is a multiple of 32, so they share a bank (32 banks, 4 bytes each), and every one of those stores is a two-way bank conflict, by construction.

The shipped code keeps the conflict on purpose. Deriving it, finding it in the report, and removing it with padding applies day 15's method to this kernel.

The honest target

The score for this day is a percentage of cuBLAS, and that number is only meaningful if both sides do the same arithmetic. The program calls cublasGemmEx with CUDA_R_32F in and out and CUBLAS_COMPUTE_32F, under CUBLAS_DEFAULT_MATH, which is plain FP32 with no TF32 tensor-core path (math modes: https://docs.nvidia.com/cuda/cublas/index.html#cublasmath-t , checked 2026-08-29). This compute type prevents newer GPUs from using a TF32 tensor-core path, so both kernels perform the same arithmetic.

Use cuBLAS on the same GPU as the denominator. The course pass mark is 50 percent of cublasGemmEx at N = 2048; the Results section records one measured comparison.

Later lessons cover three causes of the remaining gap: double buffering, asynchronous copies, and tensor cores.

Five paths, one table

Full program in matmul_registers.cu, under code/day44-matmul-2/. The harness enforces four things.

Every path computes the identical product. FP32 in, FP32 accumulate, FP32 out, checked against a double-precision CPU reference at every element of all three sizes before anything is timed, cuBLAS included.

No partial tiles, checked by static_assert. None of the kernels guards its bounds, so every case size must divide every block tile and the compiler enforces it. A bounds guard would run in the inner loop.

Production libraries often handle partial tiles with a separate kernel, but this file accepts only exact tile sizes.

The launch shape is printed before any timing. Registers per thread, spill bytes and resident blocks per SM come from cudaFuncGetAttributes and the occupancy API, which need no profiler, no counters and no root.

Events, and a warm-up per path. cuBLAS picks and caches its algorithm on the first call. Warm every path before timing it.

Step 5's inner loop is the register tiling, all of it:

        for (int p = 0; p < kBlockK; ++p) {
            for (int i = 0; i < kThreadM; ++i) {
                regM[i] = tileA[(threadRow * kThreadM + i) * kBlockK + p];
            }
            for (int j = 0; j < kThreadN; ++j) {
                regN[j] = tileB[p * kBlockN + threadCol * kThreadN + j];
            }
            for (int i = 0; i < kThreadM; ++i) {
                for (int j = 0; j < kThreadN; ++j) {
                    acc[i][j] += regM[i] * regN[j];
                }
            }
        }

Step 6's A-tile fill is where the wide load and the derived conflict live:

        // One 16-byte load pulls four consecutive k values of one row of A,
        // and the four scalar stores scatter them down one column of the
        // transposed tile. The scatter is the price of the transpose, and
        // the inner loop is what buys it back.
        for (int elem = tid; elem < (kBlockM * kBlockK) / 4;
             elem += kRegThreads) {
            const int rowA = elem / (kBlockK / 4);
            const int colA = (elem % (kBlockK / 4)) * 4;
            const float4 v = *reinterpret_cast<const float4*>(
                &a[(blockRow + rowA) * k + kIdx + colA]);
            tileA[(colA + 0) * kBlockM + rowA] = v.x;
            tileA[(colA + 1) * kBlockM + rowA] = v.y;
            tileA[(colA + 2) * kBlockM + rowA] = v.z;
            tileA[(colA + 3) * kBlockM + rowA] = v.w;
        }

And the call the whole page is scored against. cuBLAS is column major, so the program asks for C^T = B^T A^T by swapping the operands and leaving both untransposed, which lands the row-major answer in place for free:

static cublasStatus_t callCublas(cublasHandle_t handle, const float* d_a,
                                 const float* d_b, float* d_c, int size) {
    const float alpha = 1.0f;
    const float beta = 0.0f;
    return cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, size, size, size,
                        &alpha, d_b, CUDA_R_32F, size, d_a, CUDA_R_32F, size,
                        &beta, d_c, CUDA_R_32F, size, CUBLAS_COMPUTE_32F,
                        CUBLAS_GEMM_DEFAULT);
}

Results

Measured. Tesla T4, driver 595.84, CUDA 12.6 (V12.6.85), built with nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -lcublas, profiled with Nsight Compute 2024.3.2 under sudo at --clock-control base. Captured 2026-09-01 on the project's verification node; all three sizes are in the page's evidence file and the N = 2048 report ships at content/profiles/optimizing-matmul-2/.

Re-verified on 2026-09-02 with CUDA 13.0 (V13.0.88), driver 580.173.02, on the same Tesla T4. All five paths passed at all three sizes. At N = 2048 the vectorised kernel reached 67.0 percent of cuBLAS instead of 67.7 percent, which does not change any conclusion.

The fresh Nsight Systems capture is listed in evidence; the counter table below remains the CUDA 12.6 Nsight Compute capture.

The N = 2048 table, CUDA events, mean of 10 runs after 3 warm-ups:

kernel                       ms   GFLOP/s  %cuBLAS   gf/out  F/byte  smem/FMA
--------------------- --------- --------- -------- -------- ------- ---------
matmulSharedTiled        25.350     677.7     11.3    257.0    3.98    2.0000
matmulRegTile1d          10.612    1618.9     26.9     65.0   15.75    1.1250
matmulRegTile2d           6.011    2858.2     47.5     33.0   31.03    0.2500
matmulRegTile2dVec4       4.219    4071.7     67.7     33.0   31.03    0.0625
cublasGemmEx              2.858    6011.9    100.0        -       -         -
Metric, at N = 2048 step 3 4 5 6
sm__throughput.avg.pct_of_peak_sustained_elapsed 75.45 44.93 44.98 62.05
gpu__compute_memory_throughput.avg.pct_of_peak_sustained_elapsed 75.45 44.93 43.71 43.84
sm__warps_active.avg.pct_of_peak_sustained_active 99.62 49.94 47.19 47.11
launch__registers_per_thread 51 74 116 128
smsp__sass_inst_executed_op_shared_ld.sum 335,544,320 100,663,296 41,943,040 16,777,216
l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_st.sum 0 0 0 2,097,152

The five claims, against the card

  1. Partly matched. Every step beats the one before at every size, and step 3 sits at 11.3 percent, far below its 25 percent cap. The vectorised kernel landed at 67.7 percent of cublasGemmEx, 2.3 points under the 70 to 80 band this page bet on. The check's gate is the curriculum's 50 percent bar, which it clears with room; the band was a bet about this card's class and this card priced it slightly lower.
  2. Did not match. The executed shared-load instructions fall in the ratio 20 : 6 : 2.5 : 1, not 32 : 18 : 4 : 1, because ptxas widened shared loads in every kernel that offered adjacent addresses, not just step 6: an LDS.64 or LDS.128 is one executed instruction carrying two or four floats. The bytes-per-FMA column the derivation actually cares about is unchanged; what died is the assumption that instructions equal scalar loads. Day 46 reads the widths in SASS.
  3. Matched exactly. Shared store bank conflicts are zero for steps 3 to 5 and 2,097,152 for the float4 kernel, whose transposed A fill puts two lanes 512 floats apart in the same bank on every scalar store. One derivable conflict, present only where the derivation says it must be.
  4. Matched. The two 128 by 128 kernels post the lowest achieved occupancy, 47.19 and 47.11 percent against 99.62 for step 3, and the two best times. Residency is register-limited at 116 and 128 registers per thread, and it does not matter: each warp carries 64 times the work per stall.
  5. Matched. The 128 by 128 kernels score their worst %cuBLAS at N = 512: 27.1 and 32.3 there, against 34.4 and 57.5 at 1024 and 47.5 and 67.7 at 2048. A 4 by 4 grid of blocks on a 40-SM card leaves most of the machine idle no matter how good the inner loop is.

The speed-of-light pair also flipped the way the tile shapes commit it to: step 3 reads 75.45 on both throughputs, memory-shaped, while step 6 reads 62.05 compute against 43.84 memory.

Do not compare your milliseconds with these results. Compare the percentage of cublasGemmEx, the instruction ratio, and the derived bank conflict.

Also compare the speed-of-light pair as the tile size grows and the main limit changes from memory to compute.

Run it yourself

Run this on a CUDA GPU with cuBLAS. Add -lineinfo for profiling, and change -arch for your GPU if needed:

nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o matmul_registers \
    matmul_registers.cu -lcublas

The CPU reference at N = 2048 is 8.6 billion multiply-adds, which rules out an embed on its own.

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-2/: optimizing-matmul-2-512.ncu-rep, optimizing-matmul-2-1024.ncu-rep, optimizing-matmul-2-2048.ncu-rep, their .details.txt and .details.csv exports, and optimizing-matmul-2.nsys-rep. Nsight Systems needs no counters and runs anywhere. With no GPU at all, start at /setup/learn-cuda-without-a-gpu.

Exercise

The micro-tile is two constants, kThreadM and kThreadN, and the thread count falls out of them. Run the program at 4 by 4, 8 by 8 and 16 by 16, and report which shape passes the 50 percent target and what the launch table says about the ones that miss.

Time: 30 to 45 minutes. Submit: the three %cuBLAS figures at N = 2048, the three regs / local B rows, and one sentence per losing shape naming the resource it ran out of.

Check: this is a measurement day. The harness still gates correctness, all five paths against the CPU reference at every size, printing the kernel, case, row and column of the first mismatch, so a broken micro-tile cannot post a score.

The number that passes is the best shape reaching at least 50 percent of cublasGemmEx, the bar the curriculum sets; per-card bands tighter than that are unmeasured until the benchmark suite runs, so quote your percentage against prediction 1 rather than against a target. Report the two misses with their reasons either way, because the misses are the lesson.

Hint 1

Both losing shapes lose in the launch table before they lose on the clock. One line of it goes wrong for 16 by 16 and a different line for 4 by 4. What does each shape cost per thread, and what does the SM have a fixed budget of?

Hint 2

A 16 by 16 micro-tile needs 256 accumulators, and a thread on every current card can hold at most 255 registers (https://docs.nvidia.com/cuda/cuda-programming-guide/05-appendices/compute-capabilities.html , Table 31, checked 2026-08-30). Where do accumulator number 256 and its neighbours go, and which column of the launch table counts that place? For 4 by 4, count loads per FMA the way the diagram does.

Solution

The register counts predict that 8 by 8 will avoid both limits. A 16 by 16 micro-tile cannot hold its 256 accumulators in 255 registers, so the compiler spills to local memory, the local B column goes nonzero, and every FMA now waits on traffic to the memory the whole design exists to avoid. 4 by 4 fits easily but only amortises each shared load across 16 FMAs instead of 64, and its 1024-thread blocks buy occupancy the kernel no longer needs.

The best micro-tile is the largest one that does not spill. Smaller tiles issue more instructions per FMA; larger tiles use local memory.

Use the register and local-memory rows to explain the timing result.

Pitfalls

The run dies with misaligned address as soon as you touch a constant. Every reinterpret_cast<float4*> in the file is a promise that the address is a multiple of 16 bytes, and the promise is kept by the sizes, not by the types: K and N whole numbers of block tiles, micro-tile edges multiples of four. Change one and the first wide load breaks the promise. The static_asserts catch the shipped constants; they cannot see a new size you add at run time.

Your percent of cuBLAS is flattering and wrong. Two common ways. You timed cuBLAS's first call, which includes picking and caching its algorithm, so warm up every path, this program's timer does. Or on an Ampere or newer card you let the compute type drift to CUBLAS_COMPUTE_32F_FAST_TF32, at which point cuBLAS uses different arithmetic and the comparison is invalid.

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 doubled the micro-tile and the kernel got slower. The register file is not elastic: 255 registers per thread is the architectural cap, and an accumulator that does not fit becomes local memory, which uses device memory. Read local B before reading the clock. Day 17 is where that column was introduced.

Occupancy fell to half and you tried to fix it. Register tiling gives each warp more work per operand fetched, so fewer resident warps are needed to hide the same latency. A higher occupancy is not always faster; shrinking the micro-tile can increase load work. Day 45 puts a number on the occupancy these kernels gave up and what it bought.

Go deeper

Next

Day 45 measures how lower occupancy affects these kernels. Day 46 checks prediction 2 in SASS by finding the FFMA blocks and 128-bit shared loads.

Day 72 adds cp.async, double buffering, and tensor cores, then repeats the cuBLAS comparison with WMMA.