I spent longer than I should have on a warp GEMM that was wrong on exactly 75% of its output elements. The cause was one pair of constants copied from a chain that looked identical and was not. This is the write-up I wanted to find while debugging it.
Context: ROCm/rocm-libraries PR #9627, adding an FP64 16×16×4 MFMA warp GEMM to ck_tile for gfx90a, gfx942 and gfx950.
The gap
v_mfma_f64_16x16x4f64 has been exposed by the ck_tile arch layer for CDNA2, CDNA3 and CDNA4 for a long time. What was missing was the warp layer above it — the attribute implementation, the type aliases and the dispatcher specialisations that let a GEMM actually reach the instruction.
Without that chain, anyone needing a double-precision GEMM has to drop back to the legacy DeviceGemmXdl API. That is a real constraint for dense linear algebra, FEM and stencil solvers, and state-vector quantum simulation — all workloads where FP64 is not a preference but a correctness requirement.
The obvious way to build the missing chain is to copy the FP32 16×16×4 chain, which has the same tile shape, and swap the data type. That is what I did, and it is where the bug is.
What went wrong
The FP32 chain distributes the C tile across lanes one way. I assumed the FP64 instruction did the same, kept kCM0PerLane=1, kCM1PerLane=4, and got a kernel that compiled cleanly, ran, and disagreed with the CPU reference GEMM on three quarters of its outputs.
Three quarters is a useful number. It is not noise, not an accumulation-order difference, not a tolerance problem. A failure that regular is a layout failure — the values are all being computed, they are simply landing in the wrong places.
The answer is in the ISA, in the same section for both generations: MI300 (CDNA3) and MI350 (CDNA4), §7.1.3.4, V_MFMA_F64_16X16X4_F64. The FP64 output rows are packed contiguously across the lanes, rather than within each lane’s own registers the way the other MFMA instructions do. So the constants invert: kCM0PerLane=4, kCM1PerLane=1.
One more detail worth recording, because it is a genuine trap: for the f64 builtins, blgp is A/B/C negation, not lane-group select. It does not mean what it means elsewhere.
The guard is not __gfx9__
The device guard is __gfx90a__ || __gfx942__ || __gfx950__, deliberately not __gfx9__. gfx908 (CDNA1) has no FP64 matrix core at all, and a broad guard compiles a fallback that quietly returns zeros.
The CMake side had the same trap in a different shape. A guard matching gfx94|gfx95 also selects deprecated gfx940 and gfx941 and other gfx95x targets, where the header’s device guard correctly excludes the FP64 path — so the test builds the zero-returning fallback and fails for a reason that has nothing to do with the kernel. Match the header’s exact target set, and verify the negative case: I confirmed that on gfx908 the test is not registered and no FP64 device code compiles.
What the patch touches
Five files: the warp attribute implementation, the type aliases in warp_gemm.hpp, the dispatcher specialisations for fp64_t at K=4 and K=16, a unit test against ck_tile::reference_gemm, and the CMake registration scoped to the three supported targets.
Runtime-verified on gfx942 (MI300X), compiling cleanly under -Werror -Weverything. Built for gfx950; the CDNA4 f64 output layout is byte-identical to CDNA3’s — both live in §7.1.3.4, and the arch layer uses one shared definition for both — so the same result is expected there, with CI to confirm.
The two things I did not fold in
Building an example/ck_tile/03_gemm FP64 instance surfaced two pre-existing gaps in shared, precision-agnostic infrastructure that nothing had ever exercised with an 8-byte element type:
ck_tile::numeric_traits<double> does not exist. The tolerance helpers in check_err.hpp hard-require numeric_traits<T>::mant. The primary template has no mant, and the existing specialisations use uint32_t masks — too narrow for a 52-bit mantissa.
A compile-time divide-by-zero in the universal GEMM pipeline. gemm_universal_pipeline_ag_bg_cr_policy.hpp computes LdsBanksWidth / (AK1 * M0 * sizeof(ADataType)). With sizeof(ADataType)=8 on the ColumnMajor + Wave64 path, AK1 collapses and the denominator becomes zero. The shared vector-load-size logic every precision depends on does not account for 8-byte element types in this specific LDS descriptor calculation.
Neither is caused by the warp-layer change, and both belong to shared code. I filed them separately as #10381 rather than growing the original PR, so the self-contained piece can land on its own and the example can be built on a fixed base afterwards.
The takeaway
When a new precision is added to a template stack that has been precision-agnostic in theory, the assumptions that were never tested surface all at once — and they surface in shared code, not in the new path. Budget for that. The kernel was the easy half.
If you are wiring up an MFMA variant, read the ISA section for that specific instruction rather than inferring the layout from the neighbouring one. Same tile shape does not mean same lane distribution.