Coalesced Matrix-Free Finite Elements in Cell-Wise Storage

Abstract: We present a GPU-oriented formulation of continuous high-order finite elements in which the redundant, cell-wise (element-local) vector is the persistent primary representation of all field data, rather than a transient stage of matrix-free operator evaluation. We prove that, given a preconditioner whose image is continuous, the entire flexible conjugate gradient iteration can be carried out exactly on this unassembled representation: a simple primal-dual pairing identity shows that all Krylov scalars computed from local data coincide with those of the assembled solve, so inter-element communication is confined entirely to the preconditioner. The required direct stiffness summation (DSS) is then realized without indirect gather-scatter, atomics, or coloring, by a dimensionally-split cascade of one-to-one face exchanges that provably accumulates edge and vertex contributions as a byproduct of sequential axis passes; unstructured macro-block interfaces and $h$-adaptive hanging nodes are handled by disjoint topological kernels and a shadow-cell wrapper that leaves the high-throughput sweeps untouched. The cell-wise storage decouples the memory layout from the mesh topology, and we exploit this freedom to benchmark blocked layouts that trade memory coalescing against element contiguity. Numerical experiments on modern GPUs demonstrate that the resulting operator evaluation and solver outperform state-of-the-art matrix-free implementations, signifficantly exceeding throughput of existing implementations.

• The paper introduces a GPU-centric reformulation of high-order FEM by replacing global degrees of freedom with persistent, redundant cell-wise storage to enable coalesced, indirection-free memory access.
• The paper demonstrates, with rigorous proofs and numerical benchmarks, that the cell-wise approach is mathematically equivalent to classical global assembly while achieving superior throughput across polynomial degrees.
• The paper presents a deterministic cascade DSS algorithm that avoids costly synchronization operations on GPUs, paving the way for architecture-optimized, scalable FEM solvers.

Coalesced Matrix-Free Finite Elements in Cell-Wise Storage: An Expert Summary

Introduction and Motivation

The paper "Coalesced Matrix-Free Finite Elements in Cell-Wise Storage" (2607.02335) introduces a GPU-centric reformulation of continuous high-order finite element methods (FEM). The approach systematically replaces the global, non-redundant degree-of-freedom (DoF) vector typically used in matrix-free FEM with a redundant, cell-local (element-wise) vector stored persistently throughout the solver—thereby decoupling memory layout from mesh topology and enabling fully coalesced, indirection-free access patterns on modern GPUs. The work offers a mathematically rigorous framework that proves the equivalence of this representation—given appropriate preconditioners—to the classical global-assembled approach, and presents a set of algorithms for communication and summation of interface data that obviate atomics, indirect gathers/scatters, and mesh coloring. Numerical benchmarks validate that these principles translate into substantial throughput gains across polynomial degrees, especially for low to moderate $p$.

Theoretical Underpinnings

At the core is a primal-dual framework making explicit the distinction between the continuous primal space (DoFs in $\mathbb{V}_{CG}$), its redundant cell-wise image ($\mathbb{V}_{\text{cell}}$), and respective duals. Instead of treating the cell-wise representation as a transient buffer for kernel operations, this work makes it the primary persistent state of all field data. Continuity constraints are imposed solely within the preconditioner via a direct stiffness summation (DSS) operator. A key duality-pairing lemma formalizes that the dot product of an unassembled dual vector and a continuous redundant vector in $\mathbb{V}_{\text{cell}}$ yields the correct global result; thus, all Krylov subspace scalars are identical to those computed on the assembled problem.

The paper rigorously proves that, with a preconditioner mapping into the constrained (continuous) subspace, the entire flexible conjugate gradient (FCG) iteration can proceed in the cell-wise storage without explicit assembly. Only the preconditioner is required to communicate across elements, and all operator applications, vector updates, and inner products remain local and fully parallel.

DSS Algorithm and Communication Avoidance

Traditional approaches incur substantial latency from indirect gathers/scatters and atomics during the assembly of shared interface DoFs, which severely limits the attainable memory bandwidth on GPUs, especially as arithmetic throughput and memory bandwidth exhibit divergent scaling trends on contemporary hardware.

This work introduces a cascade DSS algorithm using a dimensionally-split sequence of one-to-one face exchanges (effectively $d$ axis-aligned passes for a $d$-dimensional mesh). Each pass accumulates partial sums across the corresponding interfaces, and sequential application ensures edge and vertex couplings are correctly realized without explicit handling. Boundary and non-conforming (hanging node) cases are resolved by a shadow-cell wrapper, ensuring adaptivity does not burden the high-throughput paths. For meshes with block-structured topology, the approach is further generalized using specialized face, line, and vertex kernels, accounting for arbitrary block valences and orientations.

The DSS algorithm’s deterministic access pattern and statically known offsets enable hardware-level memory coalescing, with no need for atomics or race resolution.

Figure 1: Coarse grid of the mesh used.

Numerical Performance and Storage Analysis

The storage overhead is inherent to the redundant cell-wise format: at $p=1$, storage can be $8\times$ greater than the unique DoF count, rapidly decreasing to below $1.5\times$ for $p=7$ and above. Despite higher gross memory traffic, the implementation achieves higher effective throughput due to fully coalesced memory accesses and the absence of costly indirect operations. Mass and Laplace operator evaluation rates (in GDoF/s or TFLOP/s) substantially exceed state-of-the-art (e.g., libCEED BP1), with improvements growing with polynomial degree. Structured DSS achieves $\mathbb{V}_{CG}$0 of peak memory bandwidth even as $\mathbb{V}_{CG}$1 increases, with the associated communication cost scaling favorably as surface-to-volume ratio decreases.

The design leverages the Triton language for GPU programming. Triton's tile-based execution natively suits the block-oriented storage, providing portability and performance for both NVIDIA and AMD hardware, with critical contributions upstream in support of small-tile tensor core instructions.

Contrasts and Claims

The work asserts a strong and mathematically justified claim: Flexible conjugate gradient iterations on persistent cell-wise storage, using a continuous (DSS-enforced) preconditioner, are strictly equivalent to classical global-assembled Krylov solutions on $\mathbb{V}_{CG}$2—iterate for iterate and scalar for scalar. All communication piggybacks on the preconditioner, and the persistent cell-wise representation never reconstructs the global DoF vector.

A contrastive claim is that, despite the redundancy penalty, actual throughput is uniformly higher across all practical polynomial degrees compared to classical approaches and exceeds the fastest available reference implementations on identical hardware. This is attributed to a high fraction of peak memory bandwidth utilization, enabled by the cell-wise, statically blocked, coalesced memory layout uniquely optimized for GPU performance.

Implications and Future Prospects

The cell-wise approach removes all reliance on global DoF numberings and indirect connectivity maps in the iterative solve, providing a path to truly architecture-optimized high-order FEM on modern accelerators. For production codes, this opens avenues for mixed-precision preconditioning and custom memory layouts tuned for application-specific or hardware-specific performance. The theoretical framework is generic and accommodates not only standard DSS but any preconditioner mapping into the constrained subspace (e.g., overlapping Schwarz and patch smoothers).

Practically, applications in PDE-constrained optimization, time-dependent multiphysics, and adaptive mesh refinement scenarios can exploit the decoupled layout for scalable, high-order, multinode computations with minimal architectural bottleneck. The companion multigrid work shows that even hanging node communication can be avoided, delegating all communication to the hierarchical smoother. The framework is operator-agnostic; as long as the symmetry and continuity precondition holds, these results extend directly to more general operator families.

Conclusion

This paper reframes matrix-free high-order continuous FEM by leveraging persistent, redundant cell-wise storage, providing a rigorous mathematical foundation and efficient algorithms tailored for GPU execution. By eliminating indirect memory access and confining inter-element communication to a single, highly optimized DSS operation, the method achieves superior bandwidth utilization and numerically outperforms leading alternatives. The approach generalizes naturally to more sophisticated preconditioning, mixed-precision solvers, and complex adaptive or multilevel hierarchies, and constitutes a significant advance in high-performance matrix-free FEM solver design for current and future hardware.