CPU: even-odd (symmetric) factorization for the AVX basis-apply — benchmark + RFC

[mohitt31]

Following a suggestion from @jedbrown to open an issue describing thoughts and plans — here's the RFC for adding an even-odd (symmetric) factorization path to the AVX tensor contraction kernel.

Proposal

Add an even-odd factorization code path to CeedTensorContractApply_Avx in backends/avx/ceed-avx-tensor.c. For bases whose 1D interpolation/gradient matrix has the point symmetry of GLL/GL nodes (B[i,j] = ±B[n-1-i, n-1-j]), this precomputes symmetric and antisymmetric halves and uses them to halve the FMA count in the inner contraction loop.

Evidence

I've benchmarked this in mf-kernels, a standalone set of five 1D tensor-contraction kernel variants (naive, pitfall, explicit AVX2, register-blocked AVX2, and even-odd AVX2). The isolated effect of even-odd vs register-blocked AVX2 on the bare 1D contraction:

p	blocked (GFLOP/s)	even-odd (GFLOP/s)	speedup
5	20.82	21.73	1.04×
7	21.58	25.92	1.20×
9	21.09	28.36	1.34×

(AMD EPYC Zen 3, single core, AVX2+FMA, gcc 13.3, -O3 -march=native. GFLOP/s uses the standard algorithm's FLOP count so the algorithmic saving shows as a wall-clock gain, not a smaller denominator.)

The gain grows with order because the ratio of halved FMAs to fixed overhead improves. At p=9, even-odd reaches ~28 GFLOP/s on the bare contraction.

The same repo also has a full 3D sum-factorized operator benchmark (BP1 mass + BP3 Poisson, p=1..8, structured hex mesh, ~0.5M DOFs) against MFEM native partial assembly and libCEED's /cpu/self/avx/blocked and /cpu/self/xsmm/blocked backends. This establishes the contraction kernels in a realistic operator context and confirms correctness (worst-case max-abs difference vs MFEM's own output: 3.8e-15 across all variants and orders).

Where this does NOT help — honest caveats

1. Even-odd does not beat register-blocking in the full operator. Once gather/scatter, element restriction, and the pointwise quadrature-data apply are included, the operator becomes partly memory-bound and the contraction's FMA savings wash out. In my full-operator benchmark, the register-blocked variant matches or leads even-odd at every order on both BP1 and BP3 — they only converge at p=8 BP3 (33.9 vs 33.7 M DOF/s). The CeedTensorContract is only one component of the full CeedOperator pipeline.

2. LIBXSMM leads at high order. libCEED with LIBXSMM outperforms both the mf-kernels variants and the plain AVX backend at higher orders (BP1 p=8: LIBXSMM 124.9 vs mf-kernels blocked 96.2 M DOF/s). This proposal improves the plain AVX backend's contraction kernel only, not the LIBXSMM path.

3. Scope is CeedTensorContract only. Specifically CeedTensorContract_Avx_Blocked and the dispatch in CeedTensorContractApply_Avx. Not the full CeedOperator, not restriction, not the LIBXSMM backend.

4. Portable SIMD intrinsics (AVX2/FMA via immintrin.h, matching the existing backend style), not vendor-specific. This is unrelated to the SVE backend work in SVE Backend #1984.

Proposed plan

1. This RFC — confirm the scoping makes sense.
2. WIP PR adding the even-odd path to backends/avx/ceed-avx-tensor.c:
   • Precompute the symmetric/antisymmetric halves of the basis matrix at CeedBasis setup time.
   • Gate activation on (a) detected symmetry in the 1D basis matrix and (b) a minimum order threshold so it only activates where it actually wins.
   • Existing contraction behavior is untouched for non-symmetric bases or low orders.
3. Reproducible benchmark: BP1 + BP3, p=1..8, single-core, comparing the even-odd path against the current AVX backend (/cpu/self/avx/blocked).
4. Correctness: match current backend output to ~1e-15.
5. make format + make prove-all clean.

Question

Does this scoping look right? Is there a particular order range or operator type that's most important for the AVX backend's downstream users?

[jeremylt]

That looks reasonable to me at first look.

One thing to note - I would also look at BP2 and BP4. The multi-component piece comes into play for breakpoints in ranges of orders where serial vs blocked (batched) element processing is better, and I'm curious to see if any of that changes with this work.

I'd also be curious to see whot happens to perf when we compile for single precision, though we have less robust ways to test that in the libCEED example suite.

Currently the across the board recommendation is to use LIBXSMM unless you cannot, so it would be interesting to have cases where AVX would be better.

For my understanding, would this split be applicable to LIBXSMM backend as well eventually since it specifically targets the tensor contractions, which is the main difference between the different CPU backends?

[mohitt31]

Thanks @jeremylt — really helpful pointers.

BP2 + BP4. Will add them. With ncomp=3 there's 3× the basis work per element, which should push the operator further into the compute-bound regime where the contraction's FMA savings actually matter in the full operator — as opposed to BP1/BP3 where gather/scatter dilutes the gain. I'll report where the serial/blocked/even-odd breakpoints actually land rather than assume.

Single precision. Happy to add FP32 as a throughput data point. The even-odd factorization itself is precision-agnostic (the halved FMA count applies either way), but the register tile sizes will likely need adjustment for FP32. I'll anchor correctness in FP64 and treat FP32 as a perf-only measurement, clearly flagged.

AVX vs LIBXSMM. Fully agree LIBXSMM should stay the default. The goal is to improve the plain AVX path for builds without LIBXSMM — if even-odd can close some of the gap at mid-to-high order, that's a win for the fallback. I'll include the exact crossover in the benchmark.

Even-odd on LIBXSMM. In principle yes — the decomposition is backend-agnostic. You'd feed LIBXSMM two half-size GEMMs (roughly n/2 × n/2) plus a butterfly add/sub pre/post, instead of one full-size GEMM. Whether that's a net win is empirical: at those smaller sizes LIBXSMM's per-call overhead may eat the saving. I'd like to prove it on the plain AVX path first as the clean testbed, then port to xsmm as a follow-up if the numbers support it.

Next steps: BP1–BP4 benchmark (FP64 + FP32 data point), single-core crossover analysis, then the WIP PR for the AVX path. Anything you'd want prioritized in that ordering?

[jeremylt]

For context with BP 2 and 4 - a big thing that changes is ensuring we have enough data to get full SMDI usage. Blocked backends on single component problems are great for low order because the sizes are always a multiple of 8 but around order 7 that becomes less of an issue and the serial backends start to outperform. With multi component problems that breakpoint shifts, and I'd be curious to see if anything happens with this factorization (or not).

[mohitt31]

Good point — the serial/blocked breakpoint shift is worth measuring explicitly.

The even-odd factorization halves the contraction width (roughly n/2 instead of n), which means the output tile shrinks. For blocked backends that need multiples-of-8 to fill SIMD lanes, that could push the breakpoint in either direction: the smaller tiles might actually make it harder to fill AVX registers at intermediate orders, or the reduced FMA count might compensate by finishing faster and amortizing the partial-lane penalty. With ncomp=3 the data is 3× wider at the element level, so the serial path has more to chew on before it saturates — I'd expect the serial/blocked crossover to shift upward (blocked stays competitive longer), but whether even-odd moves it further or pulls it back is genuinely unclear without measurement.

I'll instrument the benchmark to sweep serial vs blocked vs even-odd across p=1..12 for both ncomp=1 (BP1/BP3) and ncomp=3 (BP2/BP4) and report where each crossover actually lands. That should give a clean answer.

[jeremylt]

Sounds good.

As a note - I like to make sure all new contributors know that the contribution guidelines require all LLM usage to be declared for PRs submitted

[mohitt31]

Noted — I've read the contribution guidelines. Thanks for flagging.