GEMM: half/bfloat16 products accumulate in the operand type; add float-accumulating GEBP kernels

Eigen's GEMM accumulates each inner product in the result scalar type — gebp_traits fixes using AccPacket = ResPacket; (Eigen/src/Core/products/GeneralBlockPanelKernel.h:587) — so a half or bfloat16 product sums k terms in an 11-bit or 8-bit accumulator. At k = 1024 with inputs in [0, 1], where nothing cancels, the bfloat16 result is off by 22 %; the same product cast to float, multiplied, and rounded once is off by 0.17 %, a quarter of a bfloat16 ulp. The narrow kernels are also slower than float on every NEON path that widens per FMA. Two paths in the tree already accumulate in float; the rest should, and the tests should then tighten to the bound that accumulation precision earns.

Accuracy

Relative Frobenius error \|C - C_{\text{ref}}\|_F / \|C_{\text{ref}}\|_F of a 64 \times k \times 64 product, C_{\text{ref}} formed in double from the same narrow inputs. "narrow" is A * B as shipped; "f32-acc" is (A.cast<float>() * B.cast<float>()).cast<T>(), i.e. what a float-accumulating kernel delivers. 13th Gen Intel(R) Core(TM) i7-13700HX, g++ (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0 -O2, default x86-64 vectorization.

scalar inputs k narrow f32-acc ratio
bfloat16 [-1,1] 16 0.52 % 0.16 % 3.2
64 0.98 % 0.17 % 5.8
256 1.9 % 0.16 % 12
1024 3.7 % 0.17 % 22
bfloat16 [0,1] 16 0.42 % 0.18 % 2.4
64 0.76 % 0.18 % 4.4
256 2.9 % 0.18 % 16
1024 22 % 0.17 % 130
half [-1,1] 16 0.065 % 0.021 % 3.1
64 0.12 % 0.021 % 5.9
256 0.23 % 0.021 % 11
1024 0.48 % 0.021 % 23
half [0,1] 16 0.052 % 0.021 % 2.5
64 0.094 % 0.022 % 4.4
256 0.19 % 0.023 % 8.5
1024 0.77 % 0.022 % 35

The [0,1] bfloat16 row is stagnation, not rounding noise: once the partial sum passes 2^8 times the term size, a term rounds away entirely. With \varepsilon_{\text{bf16}} = 2^{-8} the deterministic bound is k\,\varepsilon = 4 at k = 1024 and even the probabilistic \sqrt{k}\,\varepsilon = 0.25; neither promises a digit. With float accumulation the error is one final rounding plus k\,\varepsilon_{\text{f32}} \approx 6 \times 10^{-5}, so the result is essentially correctly rounded for any practical k.

The test suite cannot see this. test/product_test_helpers.h bounds two evaluations of the same product against each other by \lambda\sqrt{k}\,\varepsilon_{\text{operand}}\,\||A|\,|B|\|_F — the bound accumulation in the operand precision earns — so it can only confirm that two wrong answers agree, and at bfloat16 precision that bound is of the order of the result. !3027 (merged) found the triple-product bound narrowed to half overflowing to inf from n = 160, which is what prompted this issue.

Performance

Measured 2026-08-21 with g++ 13.3 -O3 -mcpu=native on an Armv8.2-A NEON board (no native bf16): float GEMM 24.9 GF/s single-thread; bfloat16 1.93 GF/s (13x slower) even with native flags, because the generic madd converts to float and back around every FMA; half 2.65 GF/s without +fp16 (9x slower) and 49.9 GF/s with it — fast, but summing in fp16. Gating packet_traits<half>::Vectorizable off instead gives 0.27 GF/s, so there is no honest fallback: the only fix is a kernel that accumulates wide.

Where the accumulator type is fixed

  • Generic: GeneralBlockPanelKernel.h:587, using AccPacket = ResPacket;; acc() stores the accumulator into the result unchanged.
  • NEON gebp_traits<half, half, …, Architecture::NEON> (arch/NEON/GeneralBlockPanelKernel.h:186): vfmaq_n_f16 accumulates in fp16; the !EIGEN_HAS_ARM64_FP16_VECTOR_ARITHMETIC branch widens to float and narrows back per FMA — the conversion cost without the accuracy.
  • RVV10 gebp_traits<half|bfloat16, …, Architecture::RVV10> (arch/RVV10/GeneralBlockPanelKernel.h:132, 182): AccPacket = PacketXh / PacketXbf.
  • Already wide: AltiVec MMA bfloat16 (arch/AltiVec/MatrixProductMMAbfloat16.h) accumulates in __vector_quad via xvbf16ger2pp and disassembles to Packet4f; OpenBLAS sbgemm under EIGEN_USE_OPENBLAS_BFLOAT16 (products/GeneralMatrixMatrix_BLAS.h:122) returns a float result buffer that is rounded once.

Proposal

  1. Generic framework. Let AccPacket differ from ResPacket in gebp_traits — Packet<float> for half/bfloat16 operands — with a narrowing acc(). The cheapest way to feed it is convert-on-pack: gemm_pack_lhs/gemm_pack_rhs already touch every operand once, so widening there makes the kernel the float kernel with a narrowing store. Cost: the packed panels double in bytes, so evaluateProductBlockingSizesHeuristic must size k_c for float — a small loss in blocking efficiency against a 13x kernel win on NEON. Widening inside madd instead keeps the panels narrow but pays a convert per FMA, which is exactly today's slow path. The same change covers GEMV and the Tensor contraction kernels, which reuse gebp_traits.
  2. Per-backend kernels on the native widening instructions, each as a gebp_traits/gebp_kernel override like the AltiVec MMA one: NEON BFDOT/BFMMLA (bf16 → f32, Armv8.6) and FMLAL/FMLAL2 (fp16 → f32, FEAT_FHM); SVE/SME BFDOT, BFMOPA, and widening FMOPA; AVX512-BF16 vdpbf16ps and AMX tdpbf16ps/tdpfp16ps; RVV vfwmacc (Zvfh) and vfwmaccbf16 (Zvfbfwma). AVX512-FP16 has native fp16 FMA but no widening dot product; there the float kernel over vcvtph2ps-widened panels from (1) is the right choice.
  3. Tests. With accumulation in float, the expected error of a narrow-result product becomes \tfrac{1}{2}\varepsilon_{\text{narrow}}\,\|C\|_F + \lambda\sqrt{k}\,\varepsilon_{\text{f32}}\,\||A|\,|B|\|_F. verifyProduct should take the accumulation scalar as a parameter and the half/bfloat16 sweeps in product_small and product_large should tighten to it, so a kernel that silently drops back to narrow accumulation fails.
  4. API, separately. The MMA and sbgemm paths compute a float result and round it at the end; exposing that as bf16 × bf16 → float without the final rounding is what users of these types actually want, and it falls out of (1) for free. Worth a ScalarBinaryOpTraits-level design of its own rather than folding into this item.

Related: #3099 (closed) (closed; its Appendix A.1 listed these kernels and deferred the framework change), !3027 (merged).