CEED / CEED/libCEED

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

Open
#2,009 6 comments 0 reactions 0 assignees View on GitHub
Dominant language
C
Stars
265
Forks
78
Avg merge
2d 8h
Merged PRs (30d)
11

Description

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](https://github.com/mohitt31/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 #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?

Contributor guide

Open the contributing guide

Research direction

Start in backends/avx/ceed-avx-tensor.c by reading CeedTensorContract_Avx_Blocked and the dispatch in CeedTensorContractApply_Avx. Review the existing AVX backend style and the linked mf-kernels benchmark before evaluating the proposed symmetry detection and order threshold. Done means the scoped path is correct, benchmarked against the current AVX backend, and make format plus make prove-all pass.

Written by the indexing model from the issue text.

Assessment

Tech stack
c
Domain
hpc, performance
Issue type
Feature
Difficulty
5/5
Estimated time
Over a week
Activity status
Active
Clarity
Clearly specified
Newbie friendliness
38/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.