avx: even-odd decomposition for tensor contractions - #2010
Conversation
Exploit centro-symmetry of GLL-to-Gauss interpolation and gradient matrices to halve the FMA count in 1D tensor contractions for the /cpu/self/avx backend. For basis matrices where T[Q-1-q][P-1-p] == +T[q][p] (symmetric, interp) or T[Q-1-q][P-1-p] == -T[q][p] (antisymmetric, grad), split into half-size even/odd matrices t_e and t_o. Fold input, perform two half-contractions using the existing blocked/remainder/single kernels, then unfold output with symmetry-aware recombination. Half-matrices are lazily computed and cached per (t_ptr, t_mode, B, J) tuple on first apply, with an 8-entry cache. Even-odd path activates only when min(B,J) >= 4 and the matrix passes the symmetry check; otherwise falls through to the standard dispatch. Correctly handles odd B (middle column in fold), odd J (middle row in unfold), and TRANSPOSE mode. Both centro-symmetric and centro-antisymmetric cases are covered.
The operator assembly path (ceed-preconditioning.c) reuses the same t pointer with different contents across assembly calls. The even-odd cache was keyed only on (t_ptr, t_mode, B, J), so a cache hit would serve stale half-matrices when the underlying data changed. Store a copy of t in each cache entry and memcmp on every hit. If the contents changed, refresh the symmetry detection and recompute the half-matrices. This fixes t566-operator (non-symmetric multi-component mass matrix assembly) which was producing values off by a factor of 3.
jeremylt
left a comment
jeremylt
left a comment
There was a problem hiding this comment.
First comments added.
The performance comparison seems to compare opt/blocked to avx/blocked? You should compare avx/blocked on main to avx/blocked on this branch. Comparing to the opt in this branch checks how much handwritten AVX instructions help against whatever the compiler emits itself.
You should compare avx/serial on main and this branch.
let's please talk human to human for discussion in this review
Rename data structures: CacheEntry -> EvenOddEntry, cache -> entries, CacheLookup -> EvenOddLookup, CachePopulate -> EvenOddPopulate. Use descriptive variable names in fold/unfold loops (idx_lower, idx_upper, idx_folded, idx_half, idx_out_lower, idx_out_upper, idx_out_mid). Replace per-call memcmp with one-time validation. The operator assembly path (ceed-preconditioning.c) reuses the same BTD_mat pointer with different contents across (comp_in, comp_out) iterations. On the first cache hit, verify contents match the stored copy. Stable basis pointers (interp_1d, grad_1d) pass validation once and are trusted for all future calls with zero overhead. Scratch buffers that change contents are permanently invalidated and fall through to the standard path.
- Rename validated -> is_validated per naming convention - Restore wrapper functions (_4_8, _8_8) for fixed-size dispatch - Replace StandardDispatch with direct wrapper calls - Use sizeof(array) in memset for cleaner VLA zeroing
|
Some of these comments still feel LLM generated. Can you please not use any LLMs to generate any replies to my human generated questions and comments? You are permitted to use LLMs in your personal development process as long as you disclose all of the usage, but in my personal review process I need to talk to the human who is responsible for the changes to make sure we're making the right choices for the codebase |
|
@mohitt31 Do you plan to continue with this PR? If not, I will go ahead and implement it with the changes requested and credit you in the commits. |
|
Hey guys, sorry for the delay. I had my mid-sem exams. I'm still working on this. 1.Add CeedSymmetryType enum And Jeremy , I got your point about LLM replies. Sorry for that. I'l write everything myself from now on. |
|
No worries, thanks for the update! |
Move the centro-symmetry detection, half-matrix computation, and even-odd tensor contraction from the AVX backend into the interface layer so every CPU backend benefits from the optimization. - Add CeedSymmetryType enum to types.h - Add lazy-init decomposition getters in ceed-basis.c - Add CeedTensorContractApplyEvenOdd in ceed-tensor.c with correct handling of antisymmetric matrices in transpose mode - Use even-odd in ref backend for interp, collocated grad, and underintegrated grad paths - Strip even-odd code from avx backend (now uses interface-level path)
| CeedSymmetryType interp_symmetry, grad_symmetry; | ||
| const CeedScalar *interp_1d_even = NULL, *interp_1d_odd = NULL; | ||
| const CeedScalar *grad_1d_even = NULL, *grad_1d_odd = NULL; | ||
|
|
||
| CeedCallBackend(CeedBasisGetEvenOddDecompositionInterp1D(basis, &interp_symmetry, &interp_1d_even, &interp_1d_odd)); | ||
| CeedCallBackend(CeedBasisGetEvenOddDecompositionGrad1D(basis, &grad_symmetry, &grad_1d_even, &grad_1d_odd)); |
There was a problem hiding this comment.
This should be called in CeedBasisCreateTensorH1 so this decomposition is only computed once.
We'll need to modify the CeedBasis_Ref object to
typedef struct {
bool is_collocated;
CeedScalar *collo_grad_1d;
CeedSymmetryType interp_symmetry;
CeedScalar *interp_1d_even, *interp_1d_odd;
CeedSymmetryType grad_symmetry;
CeedScalar *grad_1d_even, *grad_1d_odd;
CeedSymmetryType collo_grad_symmetry;
CeedScalar *collo_grad_1d_even, *collo_grad_1d_odd;
} CeedBasis_Ref;And then CeedbasisDestoryTensor_Ref would need the update
static int CeedBasisDestroyTensor_Ref(CeedBasis basis) {
CeedBasis_Ref *impl;
CeedCallBackend(CeedBasisGetData(basis, &impl));
CeedCallBackend(CeedFree(&impl->collo_grad_1d));
CeedCallBackend(CeedFree(&impl->interp_1d_even));
CeedCallBackend(CeedFree(&impl->interp_1d_odd));
CeedCallBackend(CeedFree(&impl->grad_1d_even));
CeedCallBackend(CeedFree(&impl->grad_1d_odd));
CeedCallBackend(CeedFree(&impl->collo_grad_1d_even));
CeedCallBackend(CeedFree(&impl->collo_grad_1d_odd));
CeedCallBackend(CeedFree(&impl));
return CEED_ERROR_SUCCESS;
}There was a problem hiding this comment.
ohhh, I see what you did
You stuck this on the basis object. That can work, but then we'd still want the collo_grad_1d_even and collo_grad_1d_odd here in the backend code
|
|
||
| @return An error code: 0 - success, otherwise - failure | ||
|
|
||
| @ref Backend |
|
|
||
| @return An error code: 0 - success, otherwise - failure | ||
|
|
||
| @ref Backend |
There was a problem hiding this comment.
Same here - this isn't ever used in the backends
| CEED_SYMMETRY_SYMMETRIC = 2, | ||
| /// Centro-antisymmetric: t[j][b] = -t[J-1-j][B-1-b] | ||
| CEED_SYMMETRY_ANTISYMMETRIC = 3, | ||
| } CeedSymmetryType; |
There was a problem hiding this comment.
We'll want to add a CeedSymmetryTypes array for pretty-printing
|
Note - CI requests style fixes. Also, this change needs to be summarized in CHANGELOG.md |
The sweep only covered the top-left quadrant, so pairs formed by an upper row and a right-hand column were never compared. Matrices that are centro-symmetric in one quadrant but not the other were reported as symmetric, which silently produced wrong results in the even-odd path. Sweep the top half of the rows over all columns instead, which visits every pair at least once.
Restore the minimum dimension the AVX version had before this moved to the interface. Below 4 the fold and unfold cost more than the halved contraction saves, and the extra rounding step showed up as a tolerance failure in the multigrid doctests, where the coarse basis is 6x3.
|
Should we add a way to disable this functionality? I feel like That would make benchmarking considerably easier, which we should also do prior to merging. |
I think this is a good idea. See this commit for something similar: c7c5b1a |
3 participants
Summary
Exploit centro-symmetry of GLL-to-Gauss interpolation and gradient matrices to halve FMA count in 1D tensor contractions for the
/cpu/self/avxbackend.(t_ptr, t_mode, B, J)with content validation — no setup-time cost, correctly handles assembly paths that reuse pointers with different datamin(B,J) >= 4and matrix passes symmetry check; standard path otherwiseRFC: #2009
Benchmark results
Native x86_64 (Google Cloud Shell,
gcc -O3 -march=native), 3D tensor contractions, q=p+1, 10k iterations.ncomp=1 speedup (opt/blocked baseline / avx/blocked even-odd):
ncomp=3 speedup:
Crossover at p~5-6, reaching ~1.9x at p=12. Below p=4 the even-odd path is gated off (min dim threshold).
Test plan
/cpu/self/avx/blockedand/cpu/self/avx/serialLLM usage disclosure
Used Claude (Anthropic) for implementation assistance: drafting the even-odd fold/unfold logic, cache data structures, symmetry detection, and debugging the middle-row unfolding for antisymmetric matrices with odd J. Also used for diagnosing the t566 assembly regression (stale cache when pointer reused with different contents).