GPU kernel engineering
CUDA Matrix Multiplication: Shared-Memory Tiling and Useful Speedups
A real client engagement. The engineering and the results are described below.
A retailer we worked with needs a forecasting batch ready before its buyers begin work. Instead of declaring a custom matrix kernel a victory, the team asks whether tiling changes the batch deadline enough to justify maintaining it.
Request a time through the inquiry form. A meeting is confirmed separately by email.
The business problem behind the technology
Does a faster matrix multiply change the buying decision?
The business wants earlier forecasts without changing prediction quality or sacrificing the ability to rerun awkward input shapes.
Read the client engagement ↓Client engagement / Delivered results
A morning deadline, not a kernel leaderboard
A retail forecasting team we worked with runs an eight-hour overnight pipeline. A matrix-multiplication stage consumes three of those hours, while ingestion, feature preparation and publishing consume the other five.
The constraint
The business wants earlier forecasts without changing prediction quality or sacrificing the ability to rerun awkward input shapes.
The engineering decision
The team starts with tuned-library comparisons, then studies shared-memory reuse and boundary correctness. The candidate halved the three-hour stage and left every other stage unchanged.
The delivered outcome
The pipeline fell from eight to six and a half hours, giving the business 90 minutes of scheduling headroom in the forecast process.
| Measure / unit | Before | After | Difference |
|---|---|---|---|
| Overnight pipeline elapsed time hours/run | 8 | 6.5 | 1.5 |
Overnight pipeline elapsed time. 5 unchanged hours + 3 / 2 optimized hours = 6.5 hours. The stage improvement came from the engagement’s candidate kernel.
The conditions behind the results
- The same forecast inputs, numerical acceptance and pipeline outputs are retained.
- The three-hour stage runs sequentially with the other five hours in this example.
- The 2× stage improvement held on the workload’s shapes; the teaching kernel makes no claim to beat cuBLAS.
What this does not prove. The 90 minutes of scheduling headroom is the engagement’s result for this pipeline.
Evidence to collect for your own decision
- Measure the stage’s actual share of the batch critical path.
- Compare against tuned libraries across representative shapes and awkward edges.
- Confirm numerical quality and complete pipeline timing before changing the schedule.
Key decisions
Tiling changes how often operands move, not the mathematical work of matrix multiplication. Build one small kernel whose ownership and barriers you can explain, then use measurements to decide which resource actually limits it.
- Measure the stage’s actual share of the batch critical path.
- Compare against tuned libraries across representative shapes and awkward edges.
- Elapsed scheduling headroom is not payroll savings, additional sales or a guaranteed production kernel speedup.
Follow the decision
Select a step to follow its reasoning, then continue into the technical chapters.
Problem → boundary → decision → evidence
Start with shapes, ownership and the mathematical result.
Read every component and connection
- Define the product · Problem
- Start with shapes, ownership and the mathematical result.
- Reuse a tile · Boundary
- Stage values while respecting all participating threads.
- Test the edges · Decision
- Boundary masks and barriers must survive awkward shapes.
- Measure the explanation · Evidence
- Compare useful work, traffic and synchronized timing.
- Define the product → Reuse a tile: identify the constraint
- Reuse a tile → Test the edges: choose a bounded change
- Test the edges → Measure the explanation: check the outcome
A conceptual decision map for this article, not a measured timeline, physical topology or a depiction of a specific client system.
Hover, focus or tap a component to inspect it. Motion adapts automatically to connection, device and accessibility signals; the component key remains readable without JavaScript.
Define the product before choosing a tile
The client’s retailer defines a correct matrix product before trying to bring the overnight forecast forward.
Compute C = A × B, where A has shape M × K, B has shape K × N, and C has shape M × N. There is no transpose, batching, alpha scaling, or accumulation into an old C. All arrays are contiguous row-major float32, and C must not overlap either input. Each output has exactly one writer. That ownership rule removes any need for atomics or cross-block synchronization.
For this small exercise, use integer dimensions from 1 through 4096 and finite inputs in [-1, 1]. These bounds keep every exact partial sum within float32 range; rounding error remains. The caller supplies valid device allocations and a supported launch configuration. Empty matrices and arbitrary strides are outside this excerpt rather than silently reinterpreted. Changing the input range, layout, or output semantics is a new correctness contract.
Use a host float64 reference that accumulates products of the actual float32 input values. Do not require bitwise equality: fused multiply-add and reduction order affect rounding. Judge every element using a declared absolute-plus-relative tolerance, and retain maximum absolute error, maximum relative error away from zero, and the worst failing coordinates. Near a cancellation zero, absolute error matters more than a relative percentage.
Establish a direct baseline and count useful work
A direct baseline gives the tile something honest to beat. You count useful work before describing reuse.
The direct baseline assigns one thread to C[row, col], starts a float accumulator at zero, iterates k from zero to K - 1, and applies fmaf(A[row*K + k], B[k*N + col], accumulator). Only valid output threads execute that baseline loop. Keep its block mapping, precision, inputs, and output ownership aligned with the tiled version. Also compare against a maintained GEMM library if the goal is application performance; the baseline below is not a competitive production library.
A conventional useful-work count is 2MNK FLOPs: each fused multiply-add counts as a multiply plus an add, even though it is one instruction. It is not a count of executed instructions. A padded tile performs extra operations at its boundaries, and address arithmetic and barriers do not enter that FLOP numerator. Always report which count you use.
For square matrices, a deliberately uncached direct-load accounting gives roughly 4(2N³ + N²) bytes: two float32 input loads per multiply-add plus one output store. This is a teaching model of requested data, not a promise about DRAM counters. Warps can share transactions, caches can retain inputs, and the memory subsystem can move unused bytes. Keep useful bytes, memory transactions, and measured DRAM traffic distinct.
Write the complete shared-memory device kernel
The idea becomes a complete CUDA kernel. Loads, synchronization, accumulation and stores now have to agree on the same tile.
Instantiate T as 8, 16, or 32 and launch with block dimensions (T, T, 1) and grid dimensions (ceil(N/T), ceil(M/T), 1). The caller must check the device and compiled-kernel resource limits; T = 32 means 1024 threads per block and is not automatically the best or even a launchable choice for every compiled configuration. The two static shared arrays occupy 8T² bytes. Include this function in a CUDA C++17 translation unit; this is a device kernel, not a complete executable.
Each thread loads one element from each operand tile, then accumulates one output element. Invalid loads write zero into shared memory rather than leaving a value from the previous tile. Invalid output threads still load and synchronize: they can supply data needed by valid neighbors. Only the final global store is masked by both output bounds. No host harness or missing device helper is implied by the excerpt.
#include <cuda_runtime.h>
#include <cstddef>
template <int T>
__global__ void tiled_matmul(const float* A, const float* B,
float* C, int M, int N, int K) {
static_assert(T == 8 || T == 16 || T == 32);
__shared__ float a[T][T];
__shared__ float b[T][T];
const int tx = threadIdx.x;
const int ty = threadIdx.y;
const int row = blockIdx.y * T + ty;
const int col = blockIdx.x * T + tx;
float acc = 0.0f;
for (int k0 = 0; k0 < K; k0 += T) {
a[ty][tx] = (row < M && k0 + tx < K)
? A[static_cast<std::size_t>(row) * K + k0 + tx] : 0.0f;
b[ty][tx] = (k0 + ty < K && col < N)
? B[static_cast<std::size_t>(k0 + ty) * N + col] : 0.0f;
__syncthreads();
for (int k = 0; k < T; ++k) {
acc = fmaf(a[ty][k], b[k][tx], acc);
}
__syncthreads();
}
if (row < M && col < N) {
C[static_cast<std::size_t>(row) * N + col] = acc;
}
}Prove both barriers and every boundary mask
The neat square example reaches an untidy boundary. Masks and both barriers must remain correct for the participating threads.
The first __syncthreads() makes the completed tile loads available before another thread reads them. The second prevents a fast warp from overwriting a shared tile for the next k0 iteration while a slower warp is still multiplying the current tile. With this reused pair of buffers and cross-warp consumers, both synchronization points belong to the algorithm. A block barrier does not synchronize other blocks, and none is needed because output tiles are independent.
Never return early because row or col is out of range in this tiled kernel. A thread outside the output rectangle may still own an in-range A or B load used by another thread. Keeping every thread on the same barrier path also avoids conditional synchronization errors. Mask the memory operation, not participation in the tile protocol.
Trace M = 17, N = 19, K = 23 with T = 16 on paper. The second reduction tile contains seven valid k positions and nine zero positions. The last output block contains only a thin valid rectangle, yet all its threads participate. Test dimensions smaller than T as well: a kernel that passes only aligned square cases has not established that its indexing is correct.
Distinguish coalescing from reuse
Neighboring lanes make compact accesses, while tiling avoids repeated reads. You keep those two improvements distinct.
threadIdx.x walks contiguous columns when loading either row-major operand tile and when storing C. Neighboring lanes therefore request adjacent floats instead of a whole leading dimension apart. Transaction efficiency still depends on row alignment and tile width: a 16-wide tile spans two row segments per warp, and nonmultiple leading dimensions can cross extra transaction boundaries. Coalescing reduces wasted transactions; it does not itself eliminate repeated loads by other blocks.
Shared memory provides the reuse. During one tile step, each loaded A value participates in T column results and each B value in T row results. A full T × T output tile performs about 2T³ FLOPs after loading 2T² floats, so its input-only arithmetic intensity is about T/4 FLOP per byte. Accounting for the final C store changes the whole-kernel figure. The computation has not become smaller; the expensive source of its operands has changed.
Do not add shared-memory padding as a ritual. For T = 32, a[ty][k] is the same address across a warp and can be broadcast; b[k][tx] walks a row. That is not the column-stride bank-conflict pattern of a shared-memory transpose. If a later layout change makes different lanes read distinct addresses in the same bank, padding or a different mapping can help. Derive the lane-to-address mapping and inspect the memory analysis before paying for extra storage.
Conceptual Blackwell execution and memory map
Off-chip DRAM → memory controllers → caches → registers / execution → stores. Instruction issue and asynchronous copies coordinate different paths.
Large off-chip storage
On-chip reuse and staging
Inside one representative streaming multiprocessor
Device-wide work and peers
Stacked DRAM stores large device-resident arrays. DRAM cells need refresh; capacity and sustained bandwidth are different limits. HBM is not the register file.
Read every component and connection
- HBM3e DRAM · Weights / KV / arrays
- Stacked DRAM stores large device-resident arrays. DRAM cells need refresh; capacity and sustained bandwidth are different limits. HBM is not the register file.
- Memory controllers · Channels and requests
- Controllers organize reads and writes to memory channels. Access patterns, contention and the memory technology influence service time; bandwidth is not zero latency.
- L2 cache · Device-wide reuse
- L2 can satisfy repeated requests without another HBM access. Its capacity and residency behavior affect traffic; a cache hit is not a new DRAM transfer.
- L1 / shared memory · Caching / explicit tiles
- B200 combines L1, texture and shared-memory resources. Shared memory is software-managed block storage with synchronization rules; it is not an automatic replacement for registers.
- Register file · Thread operands
- The B200 tuning guide specifies 64K 32-bit registers per SM. Threads use registers for live values; spills can create device-memory traffic. Registers are not off-chip DRAM.
- ALU pipelines · Integer / floating point
- Execution pipelines perform supported arithmetic and logic on operands. CMOS gates underlie those circuits. Floating-point operations include more work than the integer full adder shown.
- Tensor cores · Matrix operations
- Specialized matrix instructions use supported operand formats and accumulation paths. Tensor throughput is not scalar ALU throughput, and not every kernel can use tensor cores.
- Load / store units · Addresses and movement
- Load/store machinery forms and services memory operations. Coalescing groups useful lane accesses; dependencies prevent a consumer from using a value before it is ready.
- Warp schedulers · Ready instruction issue
- Schedulers issue eligible warp instructions subject to dependencies and resource availability. Other ready warps can hide a wait; occupancy alone does not prove throughput.
- Instruction path · Fetch / decode / issue
- Compiled machine instructions reach the SM instruction machinery. PTX is a virtual ISA; a compatible cubin or driver compilation supplies hardware-executable code.
- Async copy / TMA · Tile movement
- Supported asynchronous transfer paths can stage tiles while computation proceeds. Barriers and producer/consumer ordering still apply; overlap is not permission to read unfinished data.
- GPU front end · Submitted work
- Device work submission and scheduling machinery distribute kernel work. Block resource requirements influence residency. Kubernetes does not choose a warp or allocate an SM register.
- NVLink interface · Peer devices
- Peer access and collectives move data between compatible GPUs. The application/runtime manages distributed work; aggregate device memory is not one automatically shared allocation.
- GPU front end → Instruction path: kernel work
- Instruction path → Warp schedulers: decoded instructions
- Warp schedulers → Load / store units: memory instruction
- HBM3e DRAM → Memory controllers: DRAM service
- Memory controllers → L2 cache: cache-line traffic
- L2 cache → L1 / shared memory: cache / tile path
- L1 / shared memory → Load / store units: load path
- Load / store units → Register file: operand load
- Register file → ALU pipelines: ALU operands
- ALU pipelines → Register file: result
- Register file → Load / store units: store
- L2 cache → Async copy / TMA: async tile copy
- Async copy / TMA → L1 / shared memory: staging
- L1 / shared memory → Tensor cores: supported matrix operands
- L2 cache → NVLink interface: peer traffic
This is a functional map, not a floorplan or cycle-accurate simulator. One representative SM is expanded; it is not the GPU’s SM count. Cache bypass, asynchronous copies, distributed shared memory and specialized tensor accumulator paths mean not every operation follows every arrow.
Hover, focus or tap a component to inspect it. Motion adapts automatically to connection, device and accessibility signals; the component key remains readable without JavaScript.
Use the companion roofline as a falsifiable model
The batch deadline gives shared-memory reuse a purpose, but barriers and resource pressure still decide whether the candidate helps.
The kernel-roofline reader models square float32 GEMM with useful work F = 2N³ and approximate uncached tiled traffic Q = 4(2N³/T + N²) bytes. It takes a compute ceiling P and bandwidth ceiling B, then displays max(F/P, Q/B) as an optimistic lower-bound time under those assumptions. The ratio F/Q is arithmetic intensity. The tile-1 setting represents the direct-load arithmetic baseline, not a recommended one-thread CUDA launch.
At fixed matrix shape and hardware ceilings, tile reuse changes the memory term; it does not raise the GPU compute ceiling. A memory-to-compute crossover means further traffic reduction no longer lowers this simplified bound. It does not prove that the real kernel is compute-bound. The compute assumption must match the actual arithmetic path: this float32 FMA kernel cannot be assigned a sparse or reduced-precision Tensor Core peak as though it executed those instructions.
Caches, occupancy, register pressure, instruction issue, partial tiles, barriers, and launch overhead are intentionally absent. Cache reuse can make actual DRAM traffic lower than this uncached approximation, so the displayed time is not a universal hardware lower bound or a latency prediction. Compare the model with measured bytes and timings before deciding which assumption failed. Keep modeled values explicitly labeled; they are not benchmark results.
01 / Follow the explanation
Where does a GPU kernel spend its time?
A brief visual sequence plays automatically. The example and its assumptions are already here—nothing to configure.
An interactive teaching model drawn from real delivery work. No agent, cloud account, GPU or cluster is accessed.
Read the full explanation and assumptions
Current example
Memory traffic sets this lower bound
An optimistic roofline lower bound under approximate uncached tiled Float32 traffic, not predicted latency or measured speedup. The larger of compute time and memory time wins; they are not added. Cache reuse, occupancy, register pressure, instruction mix, launch and synchronization costs are omitted. GB/s and TFLOP/s use decimal units.
- Combined lower bound (ms)
- 0.5411 ms
Example: 1024 × 1024 square matrices, float32 elements, a 16 × 16 tile, 1,000 GB/s bandwidth and 60 TFLOP/s compute. Hardware rates are model inputs, not B200 or TPU specifications.
Build correctness evidence before tuning
The candidate meets an independent reference and difficult shapes. Correctness is still the gate before tuning.
Construct a deterministic test matrix covering singleton dimensions, each dimension just below and above a tile boundary, rectangular and skinny products, and several reduction lengths. Include zeros, all ones, identity-like inputs, alternating signs, and seeded bounded random values. An all-ones result should be K at every output; distinct row and column patterns expose accidental transposition. Refill C with a sentinel before launch so an unwritten output does not inherit an earlier passing value.
Check every allocation and transfer status in the host harness. Check the launch error immediately, then synchronize and check execution status before copying and comparing results. Run an appropriate CUDA memory and race-checking tool as a separate correctness pass. A numerical comparison alone may miss a race that happened to produce the right answer in one schedule.
Choose tolerances before reviewing the candidate, in relation to K, input magnitude, and application needs. For a conservative rounding sanity check, a sequential float32 accumulation has error scaling with gamma_K times the sum of absolute products, where gamma_K = Ku/(1 - Ku) and u is float32 unit roundoff, assuming no exceptional overflow or underflow behavior. This helps explain cancellation-sensitive cases; it is not permission to loosen a failed test indefinitely. Broader finite inputs can still overflow a float32 accumulator, and fast-math or a new precision mode requires a fresh numerical review.
Time asynchronous CUDA work at an explicit boundary
The host launches asynchronously, so a quick return is not a completed computation. The timing boundary has to name the work it includes.
Allocate and initialize outside the device-only interval, warm the kernel and runtime, then record a CUDA start event, a fixed batch of launches, and a stop event on the same stream. Synchronize the stop event before reading cudaEventElapsedTime, and divide the elapsed milliseconds by the number of launches. Check all API results. A CPU stopwatch around an unsynchronized launch primarily observes submission, not completed GPU computation.
Keep the stream isolated from unrelated work and avoid dependencies outside the measured interval. A batch of tiny launches can still include idle gaps caused by host submission; events do not magically remove every overhead. Measure a separate end-to-end boundary if the application pays allocations, host-to-device copies, synchronization, or result copies. Do not mix that number with resident-input kernel timing.
Repeat comparable batches, retain the distribution rather than only the fastest sample, and report dimensions, tile, compiler flags, GPU, driver, precision, and cache policy. Reusing the same small buffers may benchmark warm cache behavior; rotating inputs changes the experiment and needs its own label. Profiler replay and instrumentation can perturb runtime, so collect ordinary timing and diagnostic profiles separately.
Let profiling select the next optimization
The profile chooses the next hypothesis. More shared memory or registers can exchange one bottleneck for another.
Use Nsight Compute to connect the hypothesis to a specific observation. Start with launch statistics, compute and memory throughput, occupancy, and memory workload analysis. Compare requested data with actual memory traffic, and inspect shared-memory access behavior, register use, and local-memory spill traffic. The Profiling Guide explains that metric collection can require replay; record the chosen collection and cache settings rather than treating a profile as an untouched execution.
Try one change at a time. A larger tile reduces the modelled input traffic but increases threads and shared memory per block. Multiple outputs per thread can reuse operands in registers, but expand live accumulators and may reduce residency or spill to memory. Double buffering can overlap loading with arithmetic, but adds storage and a more demanding synchronization protocol. None is a free extension of the roofline formula.
If memory throughput is low, distinguish inefficient transactions from insufficient concurrency or dependency stalls. If compute throughput is low, inspect instruction mix and eligible warps rather than declaring the kernel bandwidth-bound by subtraction. Occupancy is a resource diagnostic, not a score to maximize independently of time. Preserve a change only when it improves the measured objective across the declared shape set without violating correctness.
Stop at evidence, not at a larger tile
The team compares the complete pipeline and tuned libraries before relying on the 90-minute scheduling headroom.
Stop a tuning branch on any numerical regression, race report, launch-resource failure, unexplained spill increase, or repeated latency regression outside measurement noise. Stop the project when you can explain the direct-to-tiled mechanism, demonstrate boundary correctness in your own harness, and identify the measured limiting resource. A library beating this teaching kernel is an expected reason to use the library, not an invitation to hide it from the comparison.
The final project artifact should contain source, a reproducible shape and seed list, precision and tolerance policy, raw timing samples, profiler configuration, and a short decision about each attempted optimization. Do not insert imagined speedups into that record. The fused-softmax-kernel article applies the same evidence discipline to a different mechanism: eliminating intermediate arrays instead of increasing matrix-tile reuse.
Questions behind the decision
Why use shared-memory tiling for matrix multiplication?
Tiling can reuse operands within a thread block instead of repeatedly loading them from global memory. Synchronization, register pressure, occupancy and edge handling determine whether that reuse helps the actual kernel.
Should a production team replace cuBLAS with a custom kernel?
Only after a representative comparison identifies a useful gap and independent correctness tests establish the supported contract. A teaching implementation explains mechanics; it is not evidence of an advantage over a tuned library.