## Summary `EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD` decides, at compile time, whether a fixed-size product uses the coeff-based evaluator or the GEMM path. Its default has only ever been measured on one machine, and the SME work in !2873 found it badly mistuned there — in both directions. This issue is to measure it across the platforms CI already covers, and to reconsider the predicate itself, which I think is the weaker half of the problem. ## Current mechanism `Eigen/src/Core/GeneralProduct.h`: ```cpp #ifndef EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD #define EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD (2 * EIGEN_GEMM_TO_COEFFBASED_THRESHOLD) #endif static constexpr bool FixedSizeCoeffBasedProduct = ProductType == GemmProduct && Rows != Dynamic && Cols != Dynamic && Depth != Dynamic && (int(Rows) + int(Cols) + int(Depth) < EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD); ``` Provenance: `EIGEN_GEMM_TO_COEFFBASED_THRESHOLD = 20` dates from 12249849b (2017), commented "obtained on a Haswell architecture". The fixed-size threshold was added in fbc0c78a6 / !2695 (closes #1791) and tuned on the same system. Neither has been checked on another target since. ## Evidence that the value is wrong off-Haswell Apple M4, SVL=512, Apple clang 21, `Matrix<S,N,N>` products, warmed to a steady clock, best of seven. Times relative to the shipping default (lower is better — below 1.0 means the default chose the slower path): | N | 3N | float | double | |---|---|---|---| | 14 | 42 | **0.17x** | **0.60x** | | 16 | 48 | **0.31x** | 0.87x | | 18 | 54 | **0.50x** | **0.62x** | | 20 | 60 | 0.60x | 1.02x | So on this target a fixed-size 14x14 float product spends ~5.9x longer than it needs to, because the default routes it to the GEMM path. Being too *high* is just as bad: at a threshold of 80, a fixed-size 24x24 or 26x26 product is ~2x slower for both scalar types. The window that is right here is roughly 56-70; the shipping default is 40. !2873 sets 64 under `EIGEN_VECTORIZE_SME` and leaves every other target alone, because I have no data for them. I would expect other backends to want different values too — the crossover is where an unrolled scalar/packet evaluator overtakes a packed micro-kernel, which depends on register count, packet width and the kernel's blocking. I just have no evidence either way, which is the point of this issue. ## Evidence that the mechanism is weak Independent of the value, four things about the predicate look wrong: 1. **It sums the dimensions, but the cost is their product.** The coeff-based evaluator fully unrolls `Rows*Cols*Depth` multiply-adds. `2x2x100` and `34x34x34` sum to 104 and 102 — indistinguishable to a sum-based rule — while unrolling to 400 and 39304 FMAs respectively. Compile time, code size and register pressure track the product, not the sum. 2. **One number serves every scalar type.** The crossover moves with packet width and with the backend's micro-kernel block. For the *runtime* knob on SME I measured fp32 preferring ~56 and fp64 ~40 — a 1.4x spread from the scalar type alone, so one value has to compromise. 3. **The default is derived from the runtime threshold at exactly 2x.** Nothing physical couples the two crossovers, and the factor that actually measured best on SME was 1.6x (64 vs 40). Worse, the derivation is a trap: retuning `EIGEN_GEMM_TO_COEFFBASED_THRESHOLD` silently moves fixed-size dispatch as well. That is exactly what bit !2873 — raising the runtime knob to 40 pushed the fixed-size one to 80 and regressed 24x24 and 26x26 by ~2x until it was pinned. 4. **It is a global macro**, so a backend cannot express a value without users redefining it, and there is no way to make it depend on the scalar type. ## Proposed work **A. Measure.** Add a benchmark that sweeps `N` for `Matrix<S,N,N>` (plus a few non-square and thin/deep shapes) with both dispatches forced, and reports the crossover per scalar type. Run it on the platforms CI already has: linux x86-64 (SSE / AVX2 / AVX512 / generic), windows x86-64, aarch64 NEON, aarch64 SME (SVL 128/256/512/2048), arm 32-bit, ppc64le (VSX and MMA), loongarch64 LSX, riscv64 RVV. That alone would tell us whether one number can serve everyone. **B. Improve the predicate.** Candidates, roughly in order of how much I'd trust them: - gate on the unrolled work `Rows*Cols*Depth` rather than `Rows+Cols+Depth`; - additionally require the operands to plausibly fit the register file, e.g. `Rows*Depth + Depth*Cols <= c * NumberOfRegisters * PacketSize`; - route the threshold through a trait next to `gebp_traits` so a backend can override it per scalar pair, defaulting to today's value so nothing changes silently. **C. Decouple the defaults** so `EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD` is not silently redefined by anyone tuning the runtime knob, whatever value it ends up with. ## Notes for whoever picks this up **A test encodes the current value.** `test/product_small.cpp` asserts the crossover with hard-coded 13x13 and 14x14 sizes, i.e. exactly `ceil(40/3)` and one below it. Any retune breaks it. Deriving the boundary from the macro instead makes the test validate the dispatch rather than the constant: ```cpp constexpr int kGemmSide = (EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD + 2) / 3; // 3N for square typedef Eigen::Matrix<double, kGemmSide - 1, kGemmSide - 1> MatrixBelowCrossover; typedef Eigen::Matrix<double, kGemmSide, kGemmSide> MatrixAtCrossover; ``` With the current default that resolves to exactly 13 and 14, so it is a no-op until someone retunes. **Raw M4 data, in case it is useful.** !2873 measured this while adding the SME double kernel but deliberately does **not** change either threshold — one machine is not enough to justify touching a shared header, which is the whole point of this issue. Apple M4, SVL=512, Apple clang 21, warmed to a steady clock, best of seven; the fixed-size sweep is relative to the shipping default: | fixed-size threshold | float N=14 | N=16 | N=18 | N=20 | N=24 | N=26 | double N=14 | N=16 | N=18 | N=24 | N=26 | |---|---|---|---|---|---|---|---|---|---|---|---| | 56 | 0.17x | 0.31x | 0.50x | 0.90x | 0.88x | 0.89x | 0.60x | 0.87x | 0.62x | 1.01x | 1.02x | | **64** | 0.17x | 0.31x | 0.50x | 0.60x | 0.87x | 0.91x | 0.60x | 0.87x | 0.62x | 1.02x | 1.01x | | 70 | 0.17x | 0.31x | 0.50x | 0.60x | 0.87x | 0.93x | 0.60x | 0.87x | 0.62x | 1.00x | 1.01x | | 80 (what 2x gives) | 0.17x | 0.31x | 0.50x | 0.61x | **1.29x** | **1.89x** | 0.60x | 0.87x | 0.62x | **1.98x** | **1.92x** | For the runtime threshold on the same machine, geometric mean over n = 7..26 relative to the default of 20: 32 -> 0.717 float / 0.793 double; **40 -> 0.607 / 0.739**; 48 -> 0.607 / 0.820; 56 -> 0.594 / 0.838; 64 -> 0.629 / 0.893. fp32 alone prefers ~56, fp64 ~40 — a 1.4x spread from the scalar type, which is the concrete evidence for point 2 above. Related: #1791, !2695, !2873.