Day 80Module 8

in-technical-review

Capstone 4: a tensor-core GEMM

Capability

Build a matrix multiply with FP16 inputs, FP32 accumulation, and correct handling when no tile divides the matrix. Then measure it against cublasGemmEx with the same input and accumulator types.

Day 44 supplies the FP32 baseline. Days 71 to 73 supply the half-precision and tensor-core operations; day 81 is not a prerequisite because the harness owns the cuBLAS call.

Retrieve, then predict

Day 44 used a 128 by 128 block tile and kept 64 accumulators in each thread. Before reading on, answer this: why can a wider thread tile make a GEMM slower even though it reuses more data?

Check your answer

The accumulators consume registers. A wider tile can cut the number of blocks that fit on an SM, and it can spill if the compiler needs more registers than the thread can hold.

Now write four predictions:

Will a WMMA kernel that reads every fragment from global memory beat the day 44 FP32 kernel at both 1024 and 2048?
How much faster will shared-memory staging make that WMMA kernel?
Which warp tile will win: 32 by 32, 32 by 64, or 64 by 64?
What will happen when K = 617 if the K loop advances by 16 without a guarded load?

Do not open the Results section yet. Record a reason beside each prediction so you can tell whether your model or your code was wrong.

Minimal model

A fragment is the part of a small matrix held across one warp. WMMA is CUDA's C++ interface for loading fragments, multiplying them on tensor cores, and storing the result.

Mixed precision uses one type for the inputs and another for the sum. This capstone stores A and B as FP16 but keeps C in FP32 because FP16 cannot represent a finite value above 65,504.

The required WMMA shape is 16 by 16 by 16. One warp multiplies two FP16 tiles and adds the result to one FP32 accumulator tile.

Fast arithmetic still needs data reuse. A warp that loads each A and B fragment from global memory repeats loads that nearby warps could share.

Worked path: choose a block and warp tile

Start with one 128 by 128 output tile and advance through K in steps of 16. The block needs a 128 by 16 tile of A and a 16 by 128 tile of B for each step.

Those two tiles contain 4,096 FP16 values, or 8,192 bytes. They produce 2 * 128 * 128 * 16 = 524,288 floating-point operations, so staging yields 64 FLOP for each byte read into shared memory.

Choose a 32 by 64 output tile per warp. Four warp rows and two warp columns cover the block, so the block uses eight warps.

Each warp owns two 16-row fragments by four 16-column fragments. That is eight accumulator fragments, and each fragment holds eight FP32 values per lane, so each lane holds 64 accumulator values.

This choice has a reason. The warp reuses its A fragments across four output fragments and its B fragments across two, while the accumulator count stays at the register budget that worked on day 44. A CUDA thread can address at most 255 registers, so the other fragments and the address math still need room.

The kernel first stages both operand tiles, then loads fragments and updates the accumulators:

        wmma::fragment<wmma::matrix_a, kWmmaDim, kWmmaDim, kWmmaDim, __half,
                       wmma::row_major>
            aFrag[kStgFragsM];
        wmma::fragment<wmma::matrix_b, kWmmaDim, kWmmaDim, kWmmaDim, __half,
                       wmma::row_major>
            bFrag[kStgFragsN];
        for (int i = 0; i < kStgFragsM; ++i) {
            wmma::load_matrix_sync(
                aFrag[i], &tileA[(warpRow * kStgWarpM + i * kWmmaDim) * kLdA],
                kLdA);
        }
        for (int j = 0; j < kStgFragsN; ++j) {
            wmma::load_matrix_sync(
                bFrag[j], &tileB[warpCol * kStgWarpN + j * kWmmaDim], kLdB);
        }
        for (int i = 0; i < kStgFragsM; ++i) {
            for (int j = 0; j < kStgFragsN; ++j) {
                wmma::mma_sync(acc[i][j], aFrag[i], bFrag[j], acc[i][j]);
            }
        }

The first two loops load A and B from shared memory. The final nested loop reuses those fragments for every output fragment the warp owns.

Padding each shared row by 16 half values changes the bank mapping used by the fragment loads. The padding costs space but avoids the bank conflicts that would make the warp wait on shared memory.

Check the worked path

Answer these before moving on:

How many warps cover the same 128 by 128 block with a 32 by 32 warp tile?
How many accumulator values does each lane hold with a 64 by 64 warp tile?
Why can low occupancy be acceptable in the staged kernel?

Check answers

The 32 by 32 tile needs 16 warps. A 64 by 64 tile has 16 accumulator fragments, so each lane holds 128 FP32 accumulator values before address and operand registers.

The staged kernel trades resident warps for data reuse. Low occupancy is only a problem if the remaining warps cannot cover instruction or memory wait time.

Compare like with like

The score uses cuBLAS with CUDA_R_16F inputs, CUDA_R_32F output, and CUBLAS_COMPUTE_32F. Comparing this kernel with an FP32 cuBLAS call would measure two different jobs.

The harness sets CUBLAS_DEFAULT_MATH and calls cublasGemmEx for you. This is why day 81 does not gate the capstone.

The reference runner contains four paths:

Path	Arithmetic	Data path	Purpose
`matmulRegTile2dVec4`	FP32 SIMT	shared memory	day 44 bronze baseline
`matmulWmmaGlobal`	FP16 input, FP32 sum	global memory to fragments	expose repeated loads
`matmulWmmaStaged`	FP16 input, FP32 sum	padded shared tiles	required WMMA path
`cublasGemmEx`	FP16 input, FP32 sum	library-owned	matched score baseline

The reference runner checks every element against a double-precision CPU result before it times a path. Its small integer inputs are exact in FP16, so a mismatch in that runner points to an indexing or synchronization bug.

Faded practice

Practice 1: trace one K step

For a 32 by 64 warp tile, mark which A fragments each output fragment reuses. Then do the same for B and confirm that the warp issues eight matrix operations per K step.

Practice 2: remove the tile choice

Build the reference kernel three times with (kStgWarpM, kStgWarpN) set to 32 by 32, 32 by 64, and 64 by 64. Before timing, record the compiler's register count, blocks per SM, and occupancy for each build.

Do not assume that 64 by 64 spills. NVCC 12.6.2 used 96, 132, and 204 registers for the three shapes, with zero spill loads and stores in all three builds.

Practice 3: design the K tail

Trace the last iteration when K = 617. A loop that steps by 16 reaches 608, then tries to load nine valid values and seven values past the matrix edge.

Choose one repair: zero-fill invalid entries in the staged tile, or send the tail through a guarded path. Explain why zero is safe for a matrix product.

Independent gate

Implement this function without calling cuBLAS:

void solve(const __half* a, const __half* b, float* c,
           std::size_t m, std::size_t n, std::size_t k,
           cudaStream_t stream);

A is row-major, B is column-major, and C is row-major. The harness owns the matched cublasGemmEx call and the float64 Kahan reference.

Case	m, n, k	Input	Required check
`unit`	16, 16, 16	exact small integers	bitwise match
`onefrag`	16, 8, 16	exact small integers	bitwise match
`small`	128, 128, 128	symmetric uniform	tolerance match
`odd`	611, 613, 617	symmetric uniform	tolerance match, no out-of-bounds access
`skinny`	4096, 16, 4096	symmetric uniform	tolerance match
`square`	1024, 1024, 1024	normal	tolerance match
`adversarial`	256, 256, 256	edge values	tolerance match without accumulator overflow
`perf`	4096, 4096, 4096	symmetric uniform	timed after correctness

Use rtol = 1e-2 and atol = 1e-3, widened by the course's length-dependent reduction rule for the current K. The memory budget is at most 2 * K / T + 16 global reads per thread, where T is your declared tile size.

At K = 4096, the length term printed by the report is 4 * 2^-23 * 64. Keep that arithmetic in the report so the tolerance can be checked instead of guessed.

Claim the highest tier whose full gate you meet:

Tier	Requirement on compute capability 7.5	8.0, 8.9, 12.0	9.0, 10.0
pass	all eight cases correct	all eight cases correct	all eight cases correct
bronze	faster than the day 44 FP32 kernel	same	same
silver	budget met and at least 50% of cuBLAS	at least 45%	at least 30%
gold	silver plus CV below 5% and at least 70%	at least 70%	at least 45%

Submit the correctness report, the percentage of matched cuBLAS, the three-shape sweep, and a short explanation for each prediction that failed.

The full capstone takes 8 to 12 hours from a blank file. The guided sweep and tail design take 2 to 4 hours with the reference runner.

Run the reference programs

Run these commands from the code/day80-capstone-4 directory. Build the WMMA runner first:

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

Build the instruction-level stretch path:

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

Use the architecture flag for your GPU. The programs take no arguments, and most of their wall time comes from the single-threaded CPU reference.

The 2048 case asks that CPU core to perform about 8.6 billion double-precision multiply-adds. That work exceeds Compiler Explorer's 20-second cap, so its reduced 128 and 256 cases are correctness evidence only.

The stretch path uses mma.sync and ldmatrix on sm_75. It exposes the register layout that WMMA hides:

            for (int i = 0; i < kFragsM; ++i) {
                for (int j = 0; j < kFragsN; ++j) {
                    asm volatile(
                        "mma.sync.aligned.m16n8k8.row.col.f32.f16.f16.f32 "
                        "{%0,%1,%2,%3}, {%4,%5}, {%6}, {%0,%1,%2,%3};\n"
                        : "+f"(acc[i][j][0]), "+f"(acc[i][j][1]),
                          "+f"(acc[i][j][2]), "+f"(acc[i][j][3])
                        : "r"(aFrag[i][0]), "r"(aFrag[i][1]), "r"(bFrag[j]));
                }
            }

Each lane supplies two .b32 registers for A, one for B, and four FP32 accumulators. The ldmatrix instruction moves shared-memory values into the lane layout expected by mma.sync.

Neither reference program handles a ragged edge. Their 1024 and 2048 cases divide every tile, so use them to study the data path, not as proof that your independent kernel passes the eight-case gate.

Results

Compare the recorded evidence with your predictions now. These numbers describe one run, not the score your GPU must produce.

On 2026-09-02, a Tesla T4 at compute capability 7.5 with driver 580.173.02 ran both programs under CUDA 13.0 (V13.0.88). Every path matched the CPU reference at every element for N = 1024 and N = 2048.

For scale, the published T4 limits are 8.1 TFLOPS for FP32, 65 TFLOPS for mixed FP16/FP32, and 320 GB/s of memory bandwidth. Day 30 measured 8092.1 GFLOP/s on its FP32 workload; these figures are context, not a target for this capstone.

CUDA 13 event timings fill the tables below; the Nsight Compute reports remain the older CUDA 12.6 artifacts and their replay durations are not benchmark results.

path, at N = 1024	ms	GFLOP/s	%cuBLAS (f16)
`matmulRegTile2dVec4` (day 44, FP32)	1.412	1521.0	10.6
`matmulWmmaGlobal`	1.569	1368.5	9.5
`matmulWmmaStaged`	0.487	4412.7	30.7
`cublasGemmEx` (f16 in, f32 accumulate)	0.149	14393.6	100.0

path, at N = 2048	ms	GFLOP/s	%cuBLAS (f16)
`matmulRegTile2dVec4` (day 44, FP32)	9.190	1869.3	5.9
`matmulWmmaGlobal`	6.962	2467.8	7.8
`matmulWmmaStaged`	2.084	8244.4	25.9
`cublasGemmEx` (f16 in, f32 accumulate)	0.540	31839.5	100.0

The instruction-level matmulMmaStaged stretch measured 0.362 ms and 5927.5 GFLOP/s at N = 1024, then 2.190 ms and 7845.4 GFLOP/s at N = 2048.

The CUDA 12.6 N = 2048 profiler export supplies the diagnostic table below. It does not supply benchmark timings: replay changed matmulWmmaStaged from the standalone 1.907 ms to 3.32 ms.

N = 2048 exported metric	register tiled	WMMA global	WMMA staged
Compute (SM) throughput	61.98%	17.53%	21.61%
Memory throughput	43.78%	49.68%	33.64%
Achieved occupancy	47.10%	97.30%	24.90%
Scheduler cycles with no eligible warp	59.28%	92.41%	88.86%

The profiler rows come from code/day80-capstone-4/profile/tensor-core-gemm-capstone-2048.details.txt. The requested metric sm__pipe_tensor_cycles_active.avg.pct_of_peak_sustained_active is absent from all three exports, so the evidence supports no tensor-pipe percentage.

The MMA report contains two N = 1024 launches and no N = 2048 launch. Do not use it to make a claim about the larger case.

Explain the gaps

The reference run tests the first two predictions:

Global WMMA lost to the FP32 kernel at N = 1024, 1.569 versus 1.412 ms, but won at N = 2048, 6.962 versus 9.190 ms. Kernel order changed with size.
Shared staging cut WMMA time from 6.962 to 2.084 ms at N = 2048.

Staged WMMA beat the day 44 kernel but reached only 25.9% of matched cuBLAS. That earns bronze on this run and misses the 50% silver and 70% gold bars.

The earlier day 44 run put its register-tiled kernel at 67.7% of FP32 cuBLAS, whose recorded N = 2048 time was 2.858 ms. Those numbers use FP32 arithmetic, so they explain the old baseline but do not belong in the FP16 cuBLAS score.

The instruction-level path also changed order with size. It beat staged WMMA at N = 1024, 0.362 versus 0.487 ms, but lost at N = 2048, 2.190 versus 2.084 ms.

The evidence does not settle the warp-tile winner because it contains compile data for all three shapes but no matching timing sweep. It also does not run K = 617, so only the independent odd case can settle the tail prediction.

For every failed prediction, name the assumption that failed. Use matrix size, data reuse, register use, and scheduler eligibility as evidence; do not replace an explanation with the word "overhead."

Common failures

If your kernel beats cuBLAS by a large margin, check the comparison first. Pin CUBLAS_COMPUTE_32F and CUBLAS_DEFAULT_MATH, and use FP16 inputs on both sides.

If powers of two pass but 617 fails, inspect the final K step. Zero-fill the invalid part of each staged tile before any warp loads a fragment.

A misaligned address after a tile edit points to a broken load promise. Check the leading dimensions, the 16-half padding, and every float4 address.

Do not switch the accumulator type to FP16 to save registers. A 4096-term sum can overflow 65,504, and smaller sums can lose low bits long before overflow.

Sources

Day 81 shows how to make the library call that the harness used here. Your capstone result gives you the evidence needed to choose between that call and a custom kernel.