avx: even-odd decomposition for tensor contractions by mohitt31 · Pull Request #2010 · CEED/libCEED

mohitt31 · 2026-08-21T02:26:22Z

Summary

Exploit centro-symmetry of GLL-to-Gauss interpolation and gradient matrices to halve FMA count in 1D tensor contractions for the /cpu/self/avx backend.

Split symmetric (interp) / antisymmetric (grad) basis matrices into half-size even/odd components
Fold input, two half-contractions via existing blocked/remainder/single kernels, unfold output
Lazy cache (8 entries) keyed by (t_ptr, t_mode, B, J) with content validation — no setup-time cost, correctly handles assembly paths that reuse pointers with different data
Activates only when min(B,J) >= 4 and matrix passes symmetry check; standard path otherwise
Handles odd B/J, TRANSPOSE mode, both centro-symmetric and centro-antisymmetric cases

RFC: #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):

p	interp	grad	interpT	gradT
5	1.01x	1.12x	1.00x	1.06x
6	1.15x	1.15x	1.09x	1.11x
7	1.75x	1.84x	1.25x	1.58x
8	1.49x	1.58x	1.88x	1.76x
9	1.64x	1.75x	1.77x	1.61x
10	1.78x	1.60x	1.68x	1.77x
11	1.84x	1.92x	1.83x	1.84x
12	1.91x	1.89x	1.88x	1.95x

ncomp=3 speedup:

p	interp	grad	interpT	gradT
5	1.08x	1.14x	1.10x	1.12x
6	1.23x	1.21x	1.16x	1.17x
7	1.74x	1.90x	1.28x	1.60x
8	1.52x	1.61x	1.85x	1.64x
9	1.56x	1.70x	1.54x	1.62x
11	1.81x	1.88x	1.67x	1.87x
12	1.83x	1.84x	1.85x	1.77x

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

All existing basis and operator tests pass with zero regressions against both /cpu/self/avx/blocked and /cpu/self/avx/serial
Tested interp and grad in forward and transpose modes across p=2..12, dim=1..3, ncomp=1 and ncomp=3
Verified correctness for symmetric (interp_1d) and antisymmetric (grad_1d) matrices
Edge cases: odd p (middle column fold), odd q (middle row unfold), C=1 serial path
Operator assembly (t566) verified — cache validates matrix contents on every hit

LLM 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).

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

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

jeremylt · 2026-08-21T03:23:26Z

+  }
+
+  const CeedScalar tol = 1e-14;
+  int              sym = 0;


First thing - variable names should be words

Fixed, renamed variables throughout to full words.

jeremylt · 2026-08-21T03:24:53Z

+static CeedTensorContract_Avx_CacheEntry *CeedTensorContract_Avx_CacheLookup(CeedTensorContract_Avx *data, const CeedScalar *t,
+                                                                             CeedTransposeMode t_mode, CeedInt B, CeedInt J) {
+  for (CeedInt i = 0; i < data->n_cached; i++) {
+    if (data->cache[i].t_ptr == t && data->cache[i].t_mode == t_mode && data->cache[i].B == B && data->cache[i].J == J) return &data->cache[i];


What is this?

It's the even-odd entry lookup, matches a (t_ptr, t_mode, B, J) tuple against stored entries. Renamed to EvenOddLookup and added a section header.

jeremylt · 2026-08-21T03:26:19Z

+  return NULL;
+}
+
+// The t matrix contents may change between applies at the same pointer (e.g. operator assembly rewrites a scratch buffer in place), so the entry


Rewriting a scratch buffer should have no effect on application of a basis. What is the actual problem that happened so we can directly fix it?

Right, basis matrices (interp_1d, grad_1d) have stable contents so pointer-keyed lookup is fine for those.

The issue was the assembly path in ceed-preconditioning.c (line 1006). It calls CeedTensorContractApply with a scratch BTD_mat buffer that gets rewritten each (comp_in, comp_out) iteration. In t566 (p=3, q=3, dim=2, num_comp=2), the first 9x9 BTD_mat happened to be centro-symmetric, so the entry stored half-matrices for it. Next iteration same pointer, different contents, stale hit.

Replaced the per-call memcmp with one-time validation: on the first hit, verify contents still match. Stable pointers pass once and are trusted after that (zero overhead). Scratch buffers that change contents get permanently invalidated.

was this comment written by a human as I requested?

@zatkins-dev design question for you here - something about this coupling of the TensorContraction back up to its Basis so explicitly doesn't sit right but I don't see any obvious alternative

First off, honestly, I used Claude to help write my last reply. That was my mistake, and I will just use my own words from now on.

About the t566 test failure (where results were exactly 1/3 of what they should be): the issue was in ceed-preconditioning.c around line 1006.

The BTD_mat buffer is reused in a loop. My cache was only checking the memory pointer, so it accidentally kept using old cached data even when the actual data inside the buffer changed.

I fixed this by replacing the per-call data check with a one-time validation on the first cache hit.

If the data matches on that first try, it is marked as trusted (zero overhead). If it fails, it completely skips the cache. This keeps static matrices (like interp_1d) fast, while safely ignoring scratch buffers (like BTD_mat).

I'm inclined to say that it's alright, since the backend controls both. I agree that the structure is a bit odd. I need to look closer at the current relationship between the basis and tensor contraction. Is is one to one? Or one contraction per Ceed context?

I'll look in detail tomorrow.

yea, that makes sense now.... but I feel like there has to be a simpler solution.

I think at its core the complication is that we use a contraction for Basis applications and for general purpose work. Doing aBasis application alone doesn't need this safeguard.

This dual usage is mixing poorly with the TensorContraction holding Basis data in this change, I think. But I don't see an obvious design fix.

Maybe the detection of the symmetry and related work on the Basis matrices should be in the Basis object. Could have the Basis query the TensorContraction if it supports this split and if so the Basis stores the required info and calls TensorContractionApplyEvenOdd? That keeps the data in the object that actually owns it and makes this logic in the TensorContraction simpler.

Just to add a bit of context here: right now, it's exactly one tensor contraction per basis. The cache lives inside the contraction's backend data, so it naturally stays scoped to that specific basis. But I'm totally open to restructuring if you see a cleaner way to handle it!

That's actually a really clean way to separate things. If the Basis handles the symmetry check and holds the half-matrices, then the TensorContraction won't need any of this caching or validation logic at all. It would just take the pre-computed even and odd matrices directly.

I can definitely prototype this! The Basis would check the symmetry right when it's created, store the t_even and t_odd matrices, and then just call a new TensorContractionApplyEvenOdd function. Would you like me to try this approach in the next revision?

jeremylt · 2026-08-21T03:27:18Z

+  // Fold input along b
+  for (CeedInt a = 0; a < A; a++) {
+    for (CeedInt b = 0; b < B / 2; b++) {
+      const CeedInt lo = (a * B + b) * C;


real words here especially please

Done. idx_lower, idx_upper, idx_folded in fold. idx_half, idx_out_lower, idx_out_upper, idx_out_mid in unfold.

jeremylt · 2026-08-21T03:28:22Z

+  CeedScalar       *t_copy;
+  CeedScalar       *t_even;
+  CeedScalar       *t_odd;
+} CeedTensorContract_Avx_CacheEntry;


Is there a more descriptive name than 'cache'?

Renamed to EvenOddEntry with n_entries and entries[] in the container struct. Functions are EvenOddLookup, EvenOddPopulate, EvenOddValidate. Macro is CEED_AVX_EVEN_ODD_MAX_ENTRIES.

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.

jeremylt

thanks, a few more comments

jeremylt · 2026-08-21T04:03:54Z

+  CeedTransposeMode t_mode;
+  CeedInt           B, J, B_half, J_half;
+  int               symmetry;
+  bool              validated;


this would be is_validated per our convention

Fixed! Renamed it to is_validated.

jeremylt · 2026-08-21T04:07:59Z

+  }
+  return sym;
 }
-static int CeedTensorContract_Avx_Single_4_8(CeedTensorContract contract, CeedInt A, CeedInt B, CeedInt C, CeedInt J, const CeedScalar *restrict t,


This function was to help the compiler recognize a fixed size in would see frequently. Is there a particular reason it was removed in this work?

My bad on removing these, I've put them back. I originally thought the compiler would just figure out the constants through inlining, but I was wrong. The compiler actually needs to see those exact sizes right at the call site to fully unroll and vectorize the loops. The wrappers lock in those numbers (like 4,8 and 8,8) so the compiler can do its job.

jeremylt · 2026-08-21T04:09:53Z

+  CeedScalar u_even[A * B_half * C], u_odd[A * B_half * C];
+  CeedScalar w_even[A * J_half * C], w_odd[A * J_half * C];
+
+  memset(u_even, 0, sizeof(u_even));


why not just declare these to be arrays of 0s instead of using raw memset?

I really wanted to use = {0}, but since A, B_half, and C are runtime values, these are Variable-Length Arrays (VLAs). C doesn't let us initialize VLAs that way - GCC just throws a "variable-sized object may not be initialized" error. So I kept the memset, but I changed it to use sizeof(array) instead of doing the math manually. Much cleaner now!

ugh, that's right

How big can these get? Calloc/Free may be a butter fit. I've been trying to keep VLAs to small allocations and indexing casts only

Really good point. If we are dealing with a 3D problem where C gets large, these VLAs could take up way too much space on the stack. I'll swap them out and use CeedCalloc and CeedFree for the main working arrays instead. I'll only keep VLAs for really small, fixed-size stuff where we know it's safe.

jeremylt · 2026-08-21T04:11:53Z

-    // Serial C=1 Case
-    CeedTensorContract_Avx_Single_4_8(contract, A, B, C, J, t, t_mode, true, u, v);
-  } else {
-    // Blocks of 8 columns


this structure that the compiler can leverage still exists... I don't see why removing this would be helpful

Put this back for the exact same reason as the wrapper functions. StandardDispatch was passing the sizes as variables, which basically hid the constants from the compiler. Now it calls the wrappers directly.

- 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

zatkins-dev · 2026-08-21T04:28:30Z

+static int CeedTensorContract_Avx_EvenOddPopulate(CeedTensorContract_Avx *data, const CeedScalar *t, CeedTransposeMode t_mode, CeedInt B, CeedInt J,
+                                                   CeedTensorContract_Avx_EvenOddEntry **entry) {
+  *entry = NULL;
+  if (data->n_entries >= CEED_AVX_EVEN_ODD_MAX_ENTRIES) return CEED_ERROR_SUCCESS;


If we're already allocating here and cache entries are pretty light, we may as well let this cache grow arbitrarily large. We should just realloc it if we need to expand path the default size.

Good catch, that makes sense. I'll change this to a dynamically growing array. I'll start with a small default size (like 4) and just double it with CeedRealloc whenever it fills up. I'll include this in my next push.

jeremylt · 2026-08-21T06:02:50Z

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 added 2 commits August 21, 2026 04:03

jeremylt requested changes Aug 21, 2026

View reviewed changes

jeremylt added performance CPU 1-In Review labels Aug 21, 2026

jeremylt assigned jeremylt and unassigned jeremylt Aug 21, 2026

jeremylt requested changes Aug 21, 2026

View reviewed changes

avx: address second round of review feedback

1cf0c13

- 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

zatkins-dev reviewed Aug 21, 2026

View reviewed changes

Uh oh!

Conversation

mohitt31 commented Aug 21, 2026

Summary

Benchmark results

Test plan

LLM usage disclosure

Uh oh!

jeremylt left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mohitt31 Aug 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mohitt31 Aug 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mohitt31 Aug 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jeremylt Aug 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mohitt31 Aug 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mohitt31 Aug 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jeremylt left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mohitt31 Aug 21, 2026 •

edited

Loading

mohitt31 Aug 21, 2026 •

edited

Loading

mohitt31 Aug 21, 2026 •

edited

Loading

jeremylt Aug 21, 2026 •

edited

Loading

mohitt31 Aug 21, 2026 •

edited

Loading

mohitt31 Aug 21, 2026 •

edited

Loading