Tiled matrix multiplication: assign outputs and count useful work

For A with shape M by K and B with shape K by N, output C has shape M by N. One block or program instance can own a rectangular output tile and accumulate its partial sums across chunks of K. Assign each valid output exactly once unless the algorithm explicitly combines partial outputs through a sound protocol. A split-K implementation is a different reduction design, with extra combination work and numerical-order consequences. Do not add overlapping output stores merely to make the grid dimensions convenient.

When M, N, or K is not divisible by the chosen tile sizes, partial tiles require masks or a separately correct edge path. Invalid input lanes can contribute zeros to the dot product if that matches the operation, while invalid output coordinates must not be stored. The loaded layout still follows the tensor's actual strides. Bounds checks cannot be removed merely because a benchmark used square aligned dimensions. Also state whether the operation computes C=AB or a more general scaled combination involving previous C, because that changes both semantics and memory traffic.

An ideal single-read model for C=AB counts each input once and each output once. It is a lower traffic model at a named external boundary, not the total number of operand accesses made by all threads. Different tiles can reload the same input, caches can satisfy some loads, and shared-memory staging adds another transfer layer. Use an instruction count, a logical tensor-byte count, and measured device-memory bytes as separate metrics. A conventional FLOP numerator divided by inconsistent byte counters can produce a misleading intensity.

Arithmetic intensity is work divided by transferred bytes at the chosen boundary. A roofline model compares that intensity with a memory-bandwidth ceiling and a compute ceiling using matching units. It does not predict an exact kernel time: launch overhead, insufficient parallelism, instruction dependencies, conversion, and synchronization can reduce achieved throughput. A tiny matrix can be latency-bound even with favorable intensity. Compare against a tuned library under the same dtype, shape, layout, and precision contract before interpreting a custom tile as an improvement.

Larger output tiles can reuse more input values, but they require more accumulators and can reduce the number of independent blocks available for a small problem. They also change register use, shared-memory needs, edge waste, and scheduling. The best shape need not be square, especially when M and N differ greatly. A blocked algorithm's mathematical correctness is independent of whether its tile saturates a particular GPU. Separate a valid tile design from a demonstrated performance result, and record the architecture and compiler when tuning.

The roofline lab has a fixed synthetic compute ceiling of 1,000 GFLOP/s. Reset, set I=8 FLOP/byte and B=80 GB/s, predict ceiling 640 GFLOP/s, then Play or Step. Separately Reset, set I=16 and B=80, predict the compute-capped ceiling 1,000, then Play or Step. These reachable settings do not change an actual matrix tile or GPU. They illustrate why improving reuse can help only until another ceiling dominates, while the displayed ridge intensity remains 12.5 FLOP/byte for the selected bandwidth.