opt - write the operator output on first touch instead of zeroing it

[hughcars]

Purpose:

opt only registers ApplyAdd, so every CeedOperatorApply first zeroes the whole output vector. This adds an opt Apply that stores the first contribution to each output entry instead, and shares one blocked restriction between input and output fields that use the same restriction. Anything it can't handle falls back to the existing zero + ApplyAdd path. Results are bitwise identical, and avx, sve, and xsmm inherit the change.

How it works: Apply runs the same element loop as ApplyAdd; only the transpose restriction of the output changes. A node shared by several elements gets one contribution from each, and only the first may be a store, but which one comes first depends on the mesh and the element order. So the scatter keeps a bitmap with one bit per output entry, cleared at the start of each Apply: a contribution to an entry whose bit is unset stores and sets the bit, and the rest add as before. The scatter keeps the loop order of the ref transpose restriction, so every entry sums the same terms in the same order, which is why the results are bitwise identical. After the element loop, only the entries whose bit is still unset, which no element touches, get zeroed. A byte per entry instead of a bit was about 3% slower; the bitmap is 1/64 the size of the output vector.

On Graviton5 (64 ranks), BP1-6, p = 1-8 geomean against main (013a5c8) with GCC 13: opt/blocked +4.5%, sve/blocked +6.1%, xsmm/blocked +6.9%.

No cell regressed, and make prove passes on the ref, opt, sve, xsmm, and CPU gen backends. I haven't run avx or measured x86. In-place Apply, which the docs rule out, now errors on opt instead of silently returning zeros. GPU backends are unchanged: the GPU gen kernels add their outputs with atomicAdd, so no element knows it contributes first.

The second commit precomputes which element makes the first contribution to each entry. The third replaces that with a bitmap of the entries written so far, which is simpler and about 0.5% slower. I'm happy to keep either.

Update (2026-10-01):

Rebased onto main after #2061, which replaces the first commit. Against main 4958ae1 on the same Graviton5 setup: opt/blocked +3.8%, sve/blocked +5.0%, xsmm/blocked +5.4%, with no significant cell regressions. make prove passes on the same backends, also with OPENMP=1.

• I kept the precomputed masks and reverted the bitmap, which was up to 5% slower at p = 4 on a mesh numbered x-fastest. So in "How it works" above, the bitmap is replaced by masks built once at setup: walking the output restriction in transpose order marks, for each block node, the lanes whose contribution comes first, and Apply zeroes the entries no element reaches before the element loop.
• Composite operators take first touch too. Setup walks the suboperators in order with one marks array and gives each its own mask, so a DoF shared by two suboperators takes its first contribution from the first of them, and entries no suboperator touches are zeroed once. Composites with a suboperator from another backend (operators at points, /cpu/self/gen), a passive output, or overlapping component layouts keep zero + ApplyAddActive. On a two-region 3D mass composite this gives +1.7% to +5.0% on opt/blocked and xsmm/blocked for p = 1, 2, 4.
• Inside an OpenMP parallel region, Apply keeps zero + ApplyAddActive, since threads may share the output memory.
• The CeedOperatorApply docs now note that threads applying into output memory they share should zero it once and use CeedOperatorApplyAdd.
• An output that is also a passive input now errors too, as in-place Apply does.
• t512, t513, and t514 cover untouched entries, a DoF shared by suboperators of different orders, overlapping component layouts, and an empty suboperator with a passive output.

LLM/GenAI Disclosure:

This was generated in large part with opus 5.5 and reviewed with sol 6 and sol 6.1, the idea and concept of not zeroing then adding rather than setting the data was approved by me, but the implementation, testing and benchmarking were assisted.

By submitting this PR, the author certifies to its contents as described by the Developer's Certificate of Origin.

[zatkins-dev]

Just my own selfish curiosity, but do you have a plot like that for the /cpu/self/gen/blocked backend? nvm, gen uses separate restriction routines

[zatkins-dev]

This is really cool and I'm surprised how much perf it gets you! I think I lead toward the precomputed approach, constant data is good, though that may change based on my follow-up question.

How would you feel about extending this to work with composite operators? We allow overloading the ApplyComposite method, which gives you the composite operator and the input/output L-vecs. The way I see it, the touched array (your slightly slower approach) could very easily be used in the composite case to wait to zero things until the absolute end of the composite application. That matters especially if you use separate operators for different mesh subdomains, for example.

You could probably also use a precomputed approach, just keeping a mark array for the whole composite operator and separate first_touch arrays for each suboperator.

If you're not interested in doing that extension, no worries! I can throw it on my backlog and probably get to it at some point.

[jeremylt]

This becomes a problem with composite, as sometimes the suboperators have a mix of unique and shared DoFs. For example, two material regions in Ratel or mixed element topology.

I'm double checking - the big question for me is if this impl for single operators preserves the ability to select overwrite vs sum active/passive outputs independently

[zatkins-dev]

I don't think it's a problem with composite necessarily -- the default impl just zeros and calls ApplyAdd on the composite Op, so it never reaches this code path.

I think my suggestion for a possible impl would work for the overlapping DoFs -- if the composite operator is tracking whether a DoF has been touched, it could pass that info into each subsequent suboperator.

[jeremylt]

Two thoughts right now

1. How is this with CPU atomics (OpenMP)? Is there any issue?

2. The first commit is I think better handled with adding a backend interface level function CeedElemRestrictionGetBlocked which creates a blocked elem restriction of the specified block size (returning self if block size is 1 if I recall the impl guts correctly) and then holding a ref to that restriction internally so the same block rstr is recycled between calls. Then naturally we can replace that whole check of creation code with the single call and we both fix opt/gen code duplication as well as the issue you're targeting. That's big enough to justify a separate PR, I think

[jeremylt]

I don't think it's a problem with composite necessarily -- the default impl just zeros and calls ApplyAdd on the composite Op, so it never reaches this code path.

I just meant that it blocked a naive extension of this directly to composite due to that issue. Yeah, we'd need to track outputs touched at the composite operator level.

[hughcars]

Overwrite vs sum: first touch only replaces Apply for an operator whose single output is active, so passive outputs never take it. Every other operator keeps zero + ApplyAddActive, and ApplyAdd and ApplyAddActive are unchanged, so overwrite and sum for active and passive outputs behave as before.

Composite (224f00c): I kept the precomputed masks and extended them as Zach suggested, with one marks array for the composite and a mask per suboperator. A DoF shared by two material regions takes its first contribution from the first and the second adds, and entries no suboperator touches are zeroed once. Composites with suboperators at points or from /cpu/self/gen, passive outputs, or overlapping component layouts keep zero + ApplyAddActive. t513 and t514 (a0e4627) cover suboperators of different orders sharing a node, overlapping layouts, and an empty suboperator with a passive output.

OpenMP (1edf2e2): the CPU backends never open a parallel region, so this only concerns apps applying operators from their own threads, one Ceed each. ApplyAdd into shared memory doesn't go through this code. Apply into one shared output races on main too, since every thread zeroes it, so inside a parallel region opt's Apply now does what main does, and e77d462 notes in the CeedOperatorApply docs to use ApplyAdd there. In an 8-thread test, ApplyAdd into a shared array matches the serial result exactly, and per-thread Apply of single and composite operators is bitwise identical to main and to one thread.

[jeremylt]

This is the thing that is worrying me the most - this is an ElemRestriction implementation inside of the Operator object. Every time I have broken the abstraction like this I have regretted it eventually. Is there a way to let the ElemRestriction and Vector together keep track of this all? Maybe the Vector has an overwrite mask and the ElemRestriction requests and flags it?