Assembly performance optimisations - #4438
Merged
Merged
Conversation
dof_transformation_fn/dof_transformation_right_fn now return nullptr (instead of a no-op closure) when an element needs no DOF transformation, avoiding an indirect call in the hot per-cell assembly loop. This finishes guarding every remaining call site (interpolate.h, discreteoperators.h, FunctionSpace.h, Function.h, pack.h, vtk_utils.h, and the interior-facet paths in assemble_matrix_impl.h/assemble_vector_impl.h) that was still invoking the closure unconditionally and crashing with std::bad_function_call. Also fixes a mismatched guard in assemble_matrix_impl.h's interior-facet code (checked P1T but called P0), and guards the mixed-element sub_element_fns loop in FiniteElement.h, where individual heterogeneous sub-elements can return nullptr even when the composite element needs a transformation overall. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Introduces FiniteElement::with_dof_transformation_fn and
with_dof_transformation_right_fn, which invoke a continuation with
either a directly-callable, inlinable DirectDofTransform{,Right}
(for non-mixed, non-blocked elements) or the existing type-erased
std::function closure, and updates all assembly/interpolation/packing
call sites to use them.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…s paths FunctionSpace::tabulate_dof_coordinates and tabulate_lagrange_dof_coordinates are only reachable for interpolation_ident() (point-evaluation, e.g. Lagrange) elements, which never need_dof_transformations(); the DirectDofTransform branch in with_dof_transformation_fn is therefore dead code at these two call sites. Benchmarking showed the wrapping cost ~20-30% extra wall time in these hot per-cell loops for no benefit, reproducible across repeated before/after process launches. Revert just these two call sites to the plain dof_transformation_fn() form; the assembly/ interpolation call sites where the direct path is actually reachable are unaffected. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…CI)
Two real bugs, both invisible to clang locally:
1. assemble_expression_impl.h: `if constexpr (entities.rank() == 2)`
inside a nested lambda that captures `entities` by reference is
rejected by GCC ("'entities' is not a constant expression") since
forming the constant expression now goes through a runtime
reference. Use `std::remove_cvref_t<decltype(entities)>::rank()`,
which depends only on the type.
2. assemble_matrix_impl.h/assemble_vector_impl.h: `if (P0)`/`if (P1T)`
assumed the transform kernel is always contextually convertible to
bool (true for std::function and DirectDofTransform), but the
custom_kernel demo calls these impl:: functions directly with a
plain, non-nullable lambda, which has no operator bool. Add
fem::is_transform_set(), which checks truthiness only when the
callable supports it and otherwise always returns true, and use it
at all four guard sites.
Verified: reproduced both failures locally against the real headers
with Homebrew g++-16 (matching the AlmaLinux/ci-build/docs CI
failures) and icpx-style rejection of the demo (matching the oneapi
job), confirmed both are fixed, and confirmed the custom_kernel demo
still runs correctly. Full local C++ (27532 assertions) and Python
fem+io (2624 tests) suites pass.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The "Update geometry fetch" commit added compile-time block-size
dispatch to assemble_cells_matrix using
std::integral_constant<std::size_t, N>, but the function's internal
BC-zeroing loops use `for (int k = 0; k < bs0; ++k)`. Comparing an int
against an unsigned integral_constant value tripped
-Werror=sign-compare on GCC (matching assemble_vector_impl.h's
existing convention of std::integral_constant<int, N>, which doesn't
have this clash). Switch to std::integral_constant<int, N>.
Also, assemble_cells_matrix's dofmap0/dofmap1 parameters became a
constrained `auto` (DofMapPackCells) instead of a concrete
std::tuple, so the custom_kernel demo's direct calls passing bare
braced-init-lists ({dofmap.map(), 1, cells}) no longer deduce --
template arguments cannot be deduced from a braced-init-list. Wrap
them in std::tuple{...} explicitly.
Verified: cpp/test/matrix.cpp and the custom_kernel demo, which
reproduced both CI failures, now compile and the demo runs correctly.
Full local C++ suite (27532 assertions) passes.
…to garth/perf-opt
garth-wells
marked this pull request as ready for review
August 28, 2026 09:52
The raw-pointer geometry-fetch loop introduced for the matrix/vector assembly speedups compiled to three separate 8-byte scalar loads/stores per vertex instead of std::copy_n's vectorised 16+8-byte copy, causing a ~7.6% regression in vector assembly of P1 Lagrange forms (confirmed via perf-annotated production binaries, main vs branch). Restored std::copy_n for the copy itself while keeping the raw-pointer index computation. Benchmarking N1curl/RT forms (the only elements where is_direct_transform_eligible() can be true) showed DirectDofTransform gives a real win only for vector assembly of degree-1 elements and is flat-to-negative everywhere else (matrix assembly at any degree, vector assembly at degree 2). Not worth the ~200 lines in FiniteElement.h and 16 call-site conversions, so it is removed entirely; all call sites revert to the plain dof_transformation_fn() closure, keeping the independently-proven wins: nullptr-return + is_transform_set guards, and the compile-time bs dispatch in assemble_matrix_impl.h. Reverting interpolate.h/discreteoperators.h to their pre-PR structure initially reintroduced a bad_function_call crash: those functions assumed the transform closure is never null (true on main, since dof_transformation_fn() used to return a no-op closure instead of nullptr), an assumption this branch's nullptr-return change breaks. Added if (transform) guards at every previously-unconditional call site to match the convention already used correctly elsewhere. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Remote's 'Updates'/'Simplifications' commits touched the DirectDofTransform machinery (explicit lambda captures, pack_impl template simplification) that this branch has since removed entirely (see c70b518). Kept local's version on every conflicting hunk.
The merge (ba1f6c6) auto-combined remote's pack_impl signature simplification (template<int _bs,...> -> template<...>(..., auto bs, ...), no more if constexpr branching) with this branch's unmodified pack_coefficient_entity call site, since the two hunks didn't textually overlap. The result didn't compile: pack_impl<bs_c()>(...) no longer matches a signature with no non-type template parameter. Updated the call site to pass bs_c/bs as a regular argument instead of an explicit template argument, matching the new signature. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
insert_csr, insert_blocked_csr, insert_nonblocked_csr, and MatrixCSR's ghost-row unpacking each duplicated the same std::lower_bound-based column search. Factored into a shared la::impl::csr_position() helper. A hand-rolled "branchless" (conditional-move) rewrite of this search was tried first, on the theory that the column found is unpredictable from one call to the next so a branchy std::lower_bound mispredicts. It measured as a reproducible ~13% regression in scripts/bench_backend.py's csr_reasm (P1, N=100): GCC 13 did not turn the idiom into a cmov, confirmed by disassembly -- it kept a real conditional branch, structurally near-identical to libstdc++'s own lower_bound, just without libstdc++'s tuning. Reverted to a thin wrapper around std::lower_bound instead; the deduplication alone measured a genuine, reproducible +6.2% on P1/N=100 csr_reasm (no significant change on P3/N=64, where kernel cost dominates insertion cost). The failed branchless attempt is documented in a code comment so it isn't re-tried without new evidence. Drive-by: MatrixCSR::add() passed _row_ptr.size() as its debug-only row-range bound, one row looser than set()'s correct size_local()+num_ghosts(); aligned add() to match, closing a theoretical one-row out-of-bounds read for row == num_all_rows() in Developer builds. Testing: cpp/test full suite (27476 assertions) and python/test unit/la + unit/fem (2594 tests) pass unmodified. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Port the compile-time block-size dispatch and raw-pointer geometry fetch already used for cell integrals to the interior-facet and entity (exterior-facet/vertex/ridge) integral paths in assemble_matrix_impl.h and assemble_vector_impl.h, closing a gap left by the earlier cell-only optimisation. The matrix-side facet/entity kernels now take DofMapPackEntities/DofMapPackFacets-constrained dofmap arguments instead of fixed tuple types, enabling the same bs0==1/bs1==1, bs0==3/bs1==3 dispatch already used for cells. Also apply the same loop-invariant-transform-check and std::copy_n gather-loop cleanups to fem::impl::pack_impl/pack_coefficient_entity in pack.h. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
clang if-converts csr_position's std::lower_bound into a branchless cmov chain, ~35-50% slower here than GCC's branchy codegen (the cmov serializes each iteration's load behind the prior comparison, killing speculation on this short, cache-resident row). Confirmed via a standalone reproducer and disassembly; a hand-rolled loop with __builtin_expect recovers full speed, but reimplementing lower_bound over a compiler codegen quirk isn't worth it, so left as-is. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
# Conflicts: # cpp/dolfinx/fem/assemble_matrix_impl.h # cpp/dolfinx/fem/assemble_vector_impl.h
jhale
approved these changes
Aug 30, 2026
| //----------------------------------------------------------------------------- | ||
| template <std::floating_point T> | ||
| std::string FiniteElement<T>::signature() const noexcept | ||
| const std::string& FiniteElement<T>::signature() const noexcept |
Member
There was a problem hiding this comment.
High potential for dangling reference.
| /// `>= col`. Callers must check `cols[pos] == col` themselves to | ||
| /// detect a genuinely missing entry; this function does not throw. | ||
| /// | ||
| /// A single shared search used by `insert_csr`, `insert_blocked_csr`, |
Member
There was a problem hiding this comment.
Results, copy into PR text as comment.
Reverts the shared-lower_bound refactor and the add()/set() row-bound fix to their pre-PR #4438 state. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Performance optimisations for per-cell/facet/entity matrix and vector
assembly, targeting
std::functionindirection for DOF transformations andgeometry/coefficient-fetch overhead. No public API or behavioural changes.
DOF transformations: skip the closure entirely when none is needed
FiniteElement::dof_transformation_fn/dof_transformation_right_fnnowreturn
nullptrinstead of a no-op closure when an element needs no DOFtransformation, rather than allocating and invoking an empty
std::functionon every cell.invoking, via a new
fem::is_transform_set()helper (traits.h): itchecks truthiness for a nullable
std::function, and always returnstruefor a non-nullable callable such as the plain lambda thecustom_kerneldemo passes when callingimpl::assemble_*directly. Theguard is hoisted out of the per-cell/per-facet loop wherever the
transform is loop-invariant.
(non-type-erased)
DirectDofTransformpath for elements that genuinelyneed a transformation applied at assembly time (H(curl)/H(div) families
such as N1curl/RT — scalar Lagrange elements of any degree never do,
since their permutations are resolved once at dofmap-construction time,
not per cell; see
needs_dof_transformations()). Benchmarking againstN1curl/RT mass-matrix assembly (isolating just this path, holding every
other change constant) showed it only wins for vector assembly of
degree-1 elements (~4-7%) and is flat-to-negative for matrix assembly at
any degree and for vector assembly at degree 2 — not worth ~200 extra
lines in
FiniteElement.hand 16 call-site conversions, so it wasdropped. See the Benchmarks section below.
Geometry/coefficient fetch across cell, facet, and entity integrals
assemble_cells/assemble_cells_matrix,assemble_entities, andassemble_interior_facets(in bothassemble_matrix_impl.handassemble_vector_impl.h) fetch each cell's coordinate DOFs via raw-pointerindexing into the geometry dofmap and coordinate array (avoiding
md::submdspanconstruction overhead in the innermost loop), then copy thecoordinates with
std::copy_n. The same loop-invariant transform check(
is_transform_set, hoisted out of the loop) andstd::copy_ngather-loopcleanup is applied to
fem::impl::pack_impl/pack_coefficient_entityinpack.h. An earlier version of the geometry copy used a manualper-component copy loop instead of
std::copy_n, which regressed vectorassembly of P1 Lagrange forms by ~7.6%:
perf-annotated, unstrippedproduction binaries showed it compiles to three separate 8-byte scalar
loads/stores per vertex, vs.
std::copy_n's vectorised 16+8-byte copy forthe same 24 bytes — GCC 13 doesn't give a hand-written copy loop the same
treatment as
std::copy_n's specialised small-copy path. Fixed byrestoring
std::copy_nfor the copy itself while keeping the raw-pointerindex computation.
Compile-time block size for matrix assembly
impl::assemble_matrix's cell-integral, interior-facet, andexterior-facet/entity dispatch now branches on the test/trial block size
(
bs0/bs1) at compile time for the commonbs == 1andbs == 3cases(via
std::integral_constant<int, 1/3>, threaded through newDofMapPackCells/DofMapPackEntities/DofMapPackFacets-constraineddofmap parameters), falling back to the runtime
intblock sizeotherwise — mirroring the dispatch already used in
impl::assemble_vector.Other fixes found along the way
assemble_matrix_impl.h'sinterior-facet code (checked one transform but applied the other).
if constexpronentities.rank()inside a lambda that capturesentitiesby reference (rejected by GCC as not a constant expression;fixed via
remove_cvref_t<decltype(entities)>::rank()), and thesign-compare/template-deduction fallout from the
std::integral_constant-based block-size dispatch above (fixed to usestd::integral_constant<int, N>and explicitstd::tuple{...}wrappingat the
custom_kerneldemo's direct call sites).Benchmarks
Serial (1 rank), Release build, unit-cube tetrahedral mesh, median of 15
repetitions, this branch vs. current
main. Correctness cross-checked onevery run: matrix nnz/Frobenius norm identical between
mainand thisbranch in every case below.
Poisson (
-div(grad(u)) = f), scalar Lagrange:H(curl)/H(div) mass matrix (
inner(u, v)*dx) — the element familiesthat actually exercise the DOF-transformation-guard change, since Lagrange
elements never need a transformation at assembly time:
Vector-assembly deltas above are all within typical run-to-run noise for
this benchmark harness (±1-2%) except where noted; the matrix-assembly wins
(driven by the geometry-fetch fix and, for P1/P3, the new
bsdispatch)are the consistent, reproducible result of this PR. The interior-facet and
exterior-facet/entity paths (DG/discontinuous forms, Neumann/Robin
boundary integrals) receive the same geometry-fetch and
bs-dispatchtreatment but were not separately re-benchmarked; the numbers above are
representative of the cell-integral case, which dominates assembly cost
for the forms exercised here.
Testing
cpp/test): 27476 assertions pass (70test cases pass, 3 skipped as parallel-only when run serially).
fem(2384 passed) andio(198 passed; ADIOS2-dependenttests skip in this environment, unrelated to this change) test suites
pass.
custom_kerneldemo (which callsimpl::assemble_*directly with anon-nullable lambda kernel, exercising the
is_transform_setfallbackpath) builds and runs correctly, with
Assembler0(std::function-based)and
Assembler1(direct-lambda) assembly producing matching norms.