Keywords: Direct stiffness summation, matrix-free, finite elements, multigrid, GPU

AMS subject classifications: 65Y10, 65Y20, 65N55, 65N30

1 Introduction↩︎

Matrix-free evaluation of high-order finite element operators has become the method of choice for large-scale simulation on modern hardware. Instead of assembling and storing a sparse matrix, the operator action is recomputed on the fly from element-local data using sum factorization, reducing both memory footprint and memory traffic by an order of magnitude at moderate polynomial degrees [1]–[4]. This shift is driven by the hardware itself: over the past two decades the floating-point throughput of processors has grown far faster than their memory bandwidth, so that virtually all sparse and element-wise operations are limited by data movement rather than arithmetic [5]. The trend is most pronounced on GPUs, which now dominate the upper end of the TOP500 list and deliver their nominal bandwidth only under stringent conditions on the access pattern.

GPUs add two further twists. First, their memory subsystem is throughput-oriented: peak bandwidth is attained only when the threads of a warp access contiguous, aligned memory (coalescing), and any indirection-based access pattern forfeits a substantial fraction of it. Second, an ever-growing share of the silicon is devoted to matrix units (tensor cores) and to reduced-precision arithmetic, to the point that double-precision capability stagnates or is emulated through lower-precision units [6]. An algorithm aiming to maximize hardware utilization should therefore (i) arrange its data so that every load and store is coalesced, (ii) cast its inner kernels as small dense matrix products, and (iii) tolerate reduced or mixed precision inside the preconditioner [7]. These three requirements shape the design presented in this paper.

State-of-the-art matrix-free frameworks: deal.II [3], [8], [9], MFEM [10], libCEED [11], [12], Nek5000/NekRS [13], [14], and atmospheric models like NUMA [15] organize the operator action around an assembled global vector of unique degrees of freedom. Every operator application then begins by gathering element-local copies through an index map and ends by scattering the element contributions back, with atomics or mesh coloring to serialize concurrent writes. On a GPU this transition is precisely where the bandwidth is lost: the indirect, warp-divergent accesses of the gather-scatter stage are the dominant cost of the otherwise highly tuned kernels, a fact documented consistently across CEED benchmark studies [8], [16] and petascale strong-scaling analyses [14], [15]. The arithmetic of sum factorization has been polished for a decade; the remaining slack lies in how the data is stored and moved.

In this paper, we adopt an alternative perspective: we make the redundant, cell-wise (element-local) vector the persistent primary representation of all field data, and we reorganize the Krylov solver so that it operates on this representation directly. The assembled vector is never formed; continuity enters only through a direct stiffness summation (DSS) operator applied inside the preconditioner. All storage becomes block-structured with statically known offsets, so every kernel (operator application, vector updates, and DSS itself) reads and writes coalesced, indirection-free memory.

The individual ingredients of our approach (element-local cell-wise storage, sum-factorized matrix-free kernels, and direct stiffness summation) are well established. Element-local redundant vectors have been a staple of the spectral-element community since the work of Deville, Fischer, and Mund [2]; the libCEED library [11] formalizes the very decomposition we build on through its L-vector/E-vector/Q-vector hierarchy, and Nek5000/NekRS [13] store fields element-locally and invoke gather-scatter (gslib [17], [18]) whenever continuity is required. The decisive difference is that these frameworks treat the element-local vector as transient kernel input and output, while we make it the solver’s persistent state: the assembled representation is never formed at all, and all vector algebra (updates, duality pairings, and norms) is performed directly on the redundant representation.

Our first contribution is the theoretical framework that justifies this design. The key duality-pairing identity, stating that pairing an unassembled dual vector with a continuous primal vector reproduces the assembled inner product, is mathematically elementary and has long been exploited in parallel spectral-element codes and in domain-decomposition methods of FETI/BDDC type [19], [20], where inner products between one assembled and one unassembled vector are standard. What we have not found stated and proved in this form is the resulting equivalence theorem: provided the preconditioner maps into the continuous subspace, the entire flexible conjugate gradient iteration [21] executed on unassembled cell-wise data is identical, iterate by iterate, tostandard (Flexible) CG applied to the assembled global system. Assembly is thereby provably required exactly once per iteration, inside the preconditioner, and nowhere else. While standard formulations of the DSS operator rely on mapping local element-wise data to an assembled global vector [22], we bypass this requirement entirely by operating strictly on the persistent E-vector.

We note that communication is not eliminated but confined: the DSS sweep inside the smoother replaces all indirect gather-scatter traffic of the classical formulation. We also emphasize that DSS is merely the realization of this requirement adopted in this paper, not an ingredient of the theory: the equivalence theorem only demands that the preconditioner consume the residual through its assembly and return a continuous (more generally, constraint-satisfying) vector. Any preconditioner with this property qualifies, in particular overlapping patch smoothers of additive or multiplicative Schwarz type [23]–[29], whose patch-local solves produce continuous corrections by construction.

Our second contribution, and the algorithmic core of the paper, is a DSS operator that requires neither indirect addressing nor atomic operations. Existing GPU gather-scatter implementations, such as gslib on GPUs and the element restrictions of libCEED and MFEM [10], traverse index lists, typically combined with atomics or mesh coloring. We instead exploit the block-structured cell-wise layout, in the spirit of the statically addressed macro-block storage of hierarchical hybrid grids [30], [31]: a dimensionally-split cascade of three axis-wise passes, each consisting of one-to-one face exchanges with statically known offsets, provably accumulates the four-way edge and eight-way vertex sums as a byproduct of the sequential sweeps. Unstructured macro-block interfaces are handled by mutually disjoint vertex, line, and face kernels with on-the-fly integer-arithmetic orientation canonicalization. A shadow-cell wrapper that extends the same machinery to \(h\)-adaptive hanging nodes (reducing the non-conforming case exactly to the conforming one so that adaptivity never touches the high-throughput sweeps) is described in Appendix 8; it is not required by the solver, since the companion multigrid paper [32] shows that DSS across hanging nodes can be avoided altogether.

Our third contribution is an engineering study enabled by the storage paradigm itself. Because cell-wise storage decouples the memory layout from the global mesh topology, the layout becomes a free tuning parameter rather than a consequence of the assembly data structures. Blocked and interleaved element layouts have been explored before, for instance via vectorization over cells in deal.II [3], [8] and struct-of-arrays element batching in MFEM and libParanumal. However, the persistent cell-wise vector lets us benchmark the full design space, including the blocked tiling that trades memory coalescing against element contiguity, under identical solver conditions. We quantify the resulting bandwidth utilization and demonstrate, on the same A100, an assembled mass-operator throughput per unique DoF that exceeds the libCEED BP1 benchmark by up to a factor of four, the margin growing with polynomial degree.

The entire solver is implemented in Triton [33], a tile-based GPU language originating in the machine-learning community. This choice is deliberate: Triton’s programming model operates on highly efficient dense, power-of-two tiles, and takes care of coalescing, shared-memory bank conflicts, and the mapping of tile-level tensor contractions onto tensor cores, addressing concerns that would otherwise demand hand-tuned CUDA. Triton also sidesteps the CUDA toolchain entirely: its compiler translates the Python-embedded kernel source, through a sequence of intermediate representations (MLIR), down to the GPU’s native instruction set: PTX, NVIDIA’s low-level format, or its AMD counterpart; so the same kernel source runs on both NVIDIA and AMD hardware without a line of vendor-specific code. Section 5 details this implementation, including the decomposition of non-power-of-two polynomial degrees into power-of-two tiles and a small-tile tensor-core extension contributed upstream as part of this work.

The primary cost of redundancy is memory footprint, and it is concentrated at the lowest order. The overhead over the assembled vector is the asymptotic factor \(\big((p+1)/p\big)^3\), which reaches \(8\times\) at \(p=1\) but then collapses steeply with degree: to \(3.4\times\) at \(p=2\), \(2.4\times\) at \(p=3\), and toward \(1.5\times\) at \(p=7\) (Table 2), as the element interior comes to dominate the shared interface. It does not, however, translate into a runtime penalty. Because the format removes the uncoalesced gather-scatter and the destination-read that burden classical assembly, the additional bytes are streamed at the device’s coalesced peak rather than fetched through indirection, and the net effect is favourable across the entire degree range: the assembled operator already outperforms the libCEED BP1 baseline at \(p=1\) and by a widening margin thereafter (Section 6). The redundancy thus buys a larger but fully coalesced footprint in exchange for eliminating the indirect traffic that otherwise dominates the operator.

The remainder of the paper is organized as follows. Section 2 introduces the discrete setting, contrasts classical matrix-free evaluation with persistent cell-wise storage, and quantifies the memory overhead together with initial throughput measurements on the BP1 mass-operator benchmark. Section 3 develops the primal-dual framework and proves the equivalence theorem for the flexible conjugate gradient iteration. Section 4 presents the layered construction of the DSS operator: the dimensionally-split structured cascade and the unstructured interface kernels, with the shadow-cell treatment of \(h\)-adaptive non-conformity deferred to Appendix 8. Section 5 describes the implementation in the Triton language, including the mapping of sum-factorized kernels onto tensor cores and the handling of non-power-of-two polynomial degrees. Section 6 reports numerical results on an NVIDIA A100, and Section 7 concludes.

2 Continuous FEM with Cell-Wise Storage↩︎

We consider a domain \(\Omega \subset \mathbb{R}^d\) partitioned into a triangulation \(\mathcal{T}_h\) of tensor-product elements: quadrilaterals in 2D (\(d=2\)) and hexahedra in 3D (\(d=3\)). The mesh originates from a relatively coarse, conformal base grid (no hanging nodes), which fixes the topology of the domain; the computational grid is obtained from it by a hierarchical sequence of refinement cycles, either global (uniform) or local (adaptive). Local refinement inherently produces hanging nodes where cells of differing refinement levels abut. The active cells are the elements of the hierarchy without children (the finest cells in any given region), and it is on this collection that the continuous finite element space of polynomial degree \(p\) is built. The discrete problem then reads \[A u = \bar{f},\] where the global operator \(A \colon \mathbb{V}_{CG} \to \mathbb{V}_{CG}^*\) represents the discrete bilinear form acting between the continuous primal space \(\mathbb{V}_{CG}\) and its dual \(\mathbb{V}_{CG}^*\). Accordingly, \(u \in \mathbb{V}_{CG}\) and \(\bar{f} \in \mathbb{V}_{CG}^*\) are the global vectors of degrees of freedom (DoFs) and the right-hand side, respectively.

2.1 Classical matrix-free evaluation and its memory-access bottleneck↩︎

In high-performance continuous finite element frameworks, the action of \(A\) is evaluated without assembling a sparse matrix. To cleanly formalize this matrix-free workflow, we must delineate both topologies (global vs.cell-wise) and mathematical spaces (primal vs.dual). We denote the continuous global spaces as \(\mathbb{V}_{CG}\) (primal) and \(\mathbb{V}_{CG}^*\) (dual), and their redundant cell-wise counterparts as \(\mathbb{V}_{\text{cell}}\) and \(\mathbb{V}_{\text{cell}}^*\), in which every element stores its own copy of the DoFs on shared interfaces.

Transitions between these topologies in the classical matrix-free evaluation are managed by two operators:

• \(\mathcal{G} \colon \mathbb{V}_{CG} \to \mathbb{V}_{\text{cell}}\) is the global gather operator acting on primal spaces.

• \(\mathcal{G}^T \colon \mathbb{V}_{\text{cell}}^* \to \mathbb{V}_{CG}^*\) is the formal adjoint of \(\mathcal{G}\). It acts as the scatter (assembly) operator mapping local contributions back to the global topology.

With these operators established, the typical matrix-free evaluation of \(A\) is formulated as follows:

1. Gather: Extract local DoF values from the global primal vector \(u\): \(u_K = \mathcal{G}_K u\).

2. Evaluate: Compute the local unassembled dual contribution independently: \(\tilde{f}_K = A_K u_K\).

3. Scatter: Accumulate (sum) the local unassembled dual results back into the global, assembled residual vector \(\bar{f}\): \(\bar{f} = \sum_K \mathcal{G}_K^T \tilde{f}_K\).

While the arithmetic of step 2 is well understood and, with sum factorization, cheap, the cost of this pipeline on modern GPUs is dominated by steps 1 and 3. The gather and scatter are indirect memory accesses: each thread follows an element-to-global index map into a vector whose ordering is dictated by the global mesh numbering, not by the access pattern of the kernel. Consecutive threads therefore touch non-consecutive addresses, and the hardware splits each warp-wide load into many separate memory transactions.

This lack of coalescing is an important bottleneck. The nominal device bandwidth is attainable only when the threads of a warp access consecutive words that the memory controller can serve in a single transaction. ¹ An uncoalesced gather may fetch an entire cache line per useful word, discarding the bulk of the transferred data, and the scatter is worse still: concurrent accumulation into shared interface DoFs additionally requires atomic operations or mesh coloring. In practice, classical matrix-free kernels for unstructured meshes are thus limited by transaction and atomics throughput rather than by raw memory bandwidth: the headline bandwidth figure of the device is structurally out of reach.

Even if that bandwidth could be reached, it (and not the arithmetic throughput) would still be the binding constraint. The hardware landscape summarized in Table 1 explains why this paper optimizes memory traffic rather than arithmetic. For each device we list the minimum arithmetic intensity (the FLOP/byte ratio below which a kernel is memory-bound) obtained as the quotient of peak FP64 throughput and memory bandwidth. Even on the A100, the most bandwidth-friendly of the recent NVIDIA datacenter parts, a kernel must perform about five double-precision operations per transferred byte (roughly 40 per double-precision value) before the compute units matter at all; on H100-class hardware and on AMD’s CDNA accelerators the threshold lies at 10–15 FLOP/byte and beyond.

The arithmetic intensity of a sum-factorized operator is not fixed: it grows with the polynomial degree, since the work per DoF rises while the data moved per DoF does not. At low degrees the operators fall far short of the thresholds above. On the A100 and on Cartesian geometries, the kernels of this work are firmly memory-bound up to degree \(p = 2\), so the quantity that decides their performance is not the FLOP count but how fast the required bytes can be moved, i.e.the fraction of peak bandwidth actually realized.

This reframes the optimization target. The lever is not to minimize bytes moved: our cell-wise format in fact moves more data than an assembled scheme, since it stores interface DoFs redundantly. What it does instead is make every access coalesced and eliminate the read-accumulate and atomic traffic of assembly, raising the achievable bandwidth such that this gain outweighs the added redundancy. From \(p = 3\) onward this no longer holds and the kernels leave the memory-bound regime: the pipeline carries additional operations except dense matrix-matrix multiplications, so the arithmetic actually issued outstrips the useful-work count and compute, rather than bandwidth, becomes the limiter.

Although this work was initially designed to target high-order discretizations, the kernels remain memory-bound across much of the degree range; the storage strategy pays off already from degree one. Wherever bandwidth is the binding constraint, performance is governed by how the field data are laid out and accessed.

2.2 Persistent cell-wise storage↩︎

In classical implementations, the cell-wise representation \(u_K\) is never formed as a persistent global object: fragments of it are generated on-the-fly as threads read from the global vector, processed within cache, registers, or shared memory, and scattered back into the global dual vector. In this work we adopt the opposite strategy: the field data are kept persistently in the cell-wise storage format \(\mathbb{V}_{\text{cell}}\), and all mathematical operations (operator evaluations, duality pairings, vector updates, and preconditioning) are performed directly on this redundant representation. A global, non-redundant vector is never formed or traversed.

As a result, the storage layout is decoupled from the global mesh topology. Every element owns a fixed-size, self-contained set of DoFs, so these sets can be arranged in memory in whatever order suits the hardware, and every kernel addresses them through statically known offsets instead of an index map. Indirect addressing disappears from the entire solver pipeline, and with it the transaction fragmentation described above: the layout can be chosen such that warp-wide accesses are coalesced, which makes the nominal device bandwidth attainable for finite element kernels.

We remark that the framework also permits packing a small structured cluster of elements (e.g.\(2^d\) neighbouring cells) into a single logical macro-element whose interior interface DoFs are shared rather than duplicated. This reduces the storage redundancy (cf.Table 2) at the cost of slightly larger local kernels; none of the algorithms presented below depend on this packing, and we treat it as an optional compression of the same storage scheme. We do not explore it further in this work, but it is a straightforward extension of the same ideas and can be used to reduce the memory footprint of the cell-wise representation.

2.3 Initial performance measurements↩︎

To give an early indication of what the storage paradigm delivers, we measure the action of two operators on the cell-wise storage: the mass operator and the Laplace (stiffness) operator. Both act on the local DoF vector \(u_K \in \mathbb{V}_{\text{cell}}\), the cell-wise block of element \(K\) read directly from the persistent storage, and both are evaluated in sum-factorized form on each tensor-product cell, \[\begin{align} M_K u_K & = (V^T \otimes V^T \otimes V^T)\, h \,(V \otimes V \otimes V)\, u_K, \tag{1} \\ L_K u_K & = (D_x M_y M_z + M_x D_y M_z + M_x M_y D_z)\, u_K, \tag{2} \end{align}\] where \(V\) is the 1D matrix evaluating nodal values at the quadrature points, \(D\) and \(M\) are the 1D derivative and mass contractions, and the scalar \(h\) collects the cell-size scaling, mimicking a Cartesian grid and avoiding per-quadrature-point geometry data. The mass operator requires six 1D contractions, the Laplace operator seven, so the latter carries roughly three times the arithmetic; for \(p = 1, 2\) both use the even–odd decomposition and for \(p \ge 3\) they do not.

These two operators bracket the behaviour of the cell-wise format. On a naive arithmetic-intensity count one would expect the mass operator to be the worst case, since its arithmetic load is minimal and a firmly memory-bound kernel pays in full for the redundant data, as every duplicated DoF is extra traffic. However, the opposite holds: the mass contractions map poorly onto the tensor-core (MMA) tiles, so the kernel issues substantial extra arithmetic and is not, in the end, limited by memory traffic. The Laplace operator, despite its higher arithmetic intensity, tiles efficiently onto the tensor cores and sustains a high floating-point rate, so it is the operator that is memory-bound, rendering it the worst case for the cell-wise format and hence the more informative measurement.

Tables 3 and 4 report the maximum achieved throughput for the full range of polynomial degrees \(p = 1, \dots, 7\) supported by our implementation. We quote the rate both in terms of unique DoFs (the figure comparable across implementations) and cell-wise, duplicated DoFs (the figure reflecting the work actually performed), together with the attained memory throughput. The operation count per element is nearly identical to a classical matrix-free implementation, but the cell-wise variant writes its output directly instead of read-accumulating into an assembled vector. For the Laplace operator we additionally report two arithmetic rates. The theory TFLOP/s is the useful work: it multiplies the cell-wise DoF rate by the floating-point operations the sum-factorized operator of 2 must perform per DoF. The achieved TFLOP/s instead counts the operations the hardware actually issues; on the tensor-core path used for \(p \ge 3\) the contraction tiles are padded, so the hardware performs more arithmetic than is necessary and the achieved rate exceeds the theoretical one. For \(p = 1, 2\) no tensor cores and hence no padding are involved, so the two rates coincide.

The tables bear this out: at the upper degrees the Laplace operator delivers the higher throughput per unique DoF of the two (cf.Tables 3 and 4), its three-fold arithmetic absorbed almost for free by the better-filled tensor-core tiles. The reference row of Table 4 lists the achieved Laplace rates of Cui et al. [28] (reported for \(p = 2, \dots, 8\) and aligned here to the matching degree); our useful-work (theory) rate exceeds theirs at every common degree. At \(p = 2\) this margin is moreover understated: the reference figures do not employ the even–odd factorization, so they count arithmetic our formulation avoids, and our advantage in actual throughput (GDoF/s) is correspondingly larger than the TFLOP/s rows alone convey. A full benchmark campaign, including layout comparisons and the direct stiffness summation (DSS) kernels, follows in Section 6.

3 Krylov Subspace Methods with explicit Primal-Dual pairing↩︎

In a cell-wise (DG-like) storage scheme for continuous finite elements, the global degrees of freedom are duplicated across shared element interfaces and hanging nodes. This duplication might appear to be a liability: every operation that touches shared DoFs seems to require synchronization of the redundant copies. However, the opposite holds, provided one maintains a mathematical distinction between two kinds of objects living in the redundant storage: continuous primal fields, whose duplicated entries agree by construction, and broken dual vectors of local cell-wise integrals, which are meaningful without any summation. Krylov subspace methods, when written so that every inner product pairs one object of each kind, never need to convert between the two representations. We exploit this to reformulate the (flexible) conjugate gradient method so that the operator application, all vector updates, and all duality pairings are evaluated without any inter-element communication, and we prove that the resulting iteration is identical, iterate by iterate, to the classical solver applied to the assembled system. The only component that must bridge the two representations is the preconditioner, which we address in Section 3.3.

3.1 Constraint Operator and Duality Pairing Identity↩︎

Let \(\mathbb{V}_{\text{cell}}\) denote the high-dimensional broken space encompassing all local degrees of freedom stored independently on each cell. We define the continuous subspace as the image of the gather operator \(\mathcal{G}(\mathbb{V}_{CG}) \subset \mathbb{V}_{\text{cell}}\), representing continuous fields.

We introduce a constraint matrix \(C \colon \mathbb{V}_{\text{cell}} \to \mathbb{V}_{\text{cell}}\) that projects any broken vector onto the continuous subspace. The exact formulation of \(C\) can vary (e.g., averaging values on shared faces), but its defining mathematical property is that it is idempotent: \(C = C^2\), with its range being exactly \(\mathcal{G}(\mathbb{V}_{CG})\). Consequently, any vector \(u \in \mathbb{V}_{\text{cell}}\) is a valid continuous primal vector if and only if it is a fixed point of the constraint operator: \[\label{eq:constraint} u = C u.\tag{3}\] Because \(C\) can represent any idempotent projection, Eq. 3 handles complex spatial constraints, ensuring that hanging nodes are slaved to their coarse neighbors.

Conversely, the dual space \(\mathbb{V}_{\text{cell}}^*\) is the space of linear functionals (namely, the raw, unassembled local residuals and fluxes. Dual vectors, denoted \(\tilde{f} \in \mathbb{V}_{\text{cell}}^*\), represent local integrals and do not require inter-element summation to be valid. The true assembled global residual representation is recovered via the transpose of the constraint matrix, \(C^T \tilde{f}\), which physically represents the accumulation across interfaces.

The advantage of cell-wise storage is evident when evaluating the standard Euclidean duality pairing \(\langle \cdot, \cdot \rangle\) between an unassembled dual vector and a continuous primal vector.

Note that the fixed-point property of continuous fields can be written as \(C \mathcal{G} = \mathcal{G}\), since \(\text{Im}(\mathcal{G}) = \text{Im}(C)\) and \(C\) acts as the identity on its range. An immediate consequence of Lemma 1 is that the pairing is also independent of the particular choice of \(C\) and recovers the pairing on the minimal global space: for \(u = \mathcal{G} u^{\text{global}}\) with \(u^{\text{global}} \in \mathbb{V}_{CG}\), \[\label{eq:pairing95global} \langle \tilde{f}, u \rangle = \langle \tilde{f}, \mathcal{G} u^{\text{global}} \rangle = \langle \mathcal{G}^T \tilde{f}, u^{\text{global}} \rangle,\tag{4}\] i.e.the local pairing of an unassembled dual vector with a continuous primal vector equals the global pairing of the assembled functional \(\mathcal{G}^T \tilde{f}\) with the underlying global coefficient vector.

This simple identity allows us to compute the global duality pairing using only local data. Given the primal vector \(u\) satisfying the constraints, pairing it with the raw cell-wise residuals is enough to recover the global result. This allows us to calculate accurate global norms and duality pairings without ever needing to perform an explicit assembly.

With this duality established, we can formally state the governing equation of the continuous finite element problem. The true continuous problem seeks a primal vector \(u \in \mathbb{V}_{\text{cell}}\) satisfying the continuity constraint \(u = C u\) such that the globally assembled residual vanishes for any continuous test function. Let \(\tilde{b} \in \mathbb{V}_{\text{cell}}^*\) be the unassembled right-hand side. The globally assembled linear system is formally expressed as: \[\label{eq:constrained95system} C^T A C u = C^T \tilde{b},\tag{5}\] where \(C^T A C \colon \mathbb{V}_{\text{cell}} \to \mathbb{V}_{\text{cell}}^*\) represents the complete system operator embedded within the redundant space. This constrained formulation is equivalent to the classical minimal global system \(A_{CG} u^{\text{global}} = b_{CG}\) posed over the continuous space \(\mathbb{V}_{CG}\), where \(A_{CG} = \mathcal{G}^T A \mathcal{G}\) is the explicitly assembled operator, \(b_{CG} = \mathcal{G}^T \tilde{b}\), and the corresponding solutions map identically via \(u = \mathcal{G} u^{\text{global}}\).

To deploy this identity within a preconditioned Krylov solver, we must distinguish how operators map between the primal and dual spaces. Central to the problem definition is the local operator \(A\), which maps a continuous primal vector to an unassembled dual residual. In a preconditioned solver, an additional ingredient is the preconditioner \(P\), which attempts to invert \(A\) approximately by mapping a dual residual back to a primal correction. Concretely:

• The Local Operator (\(A\)): The application of the local element matrices maps a continuous primal search direction \(p\) to an unassembled dual flux \(\tilde{q}\). No inter-element communication is performed. \(A \colon \mathbb{V}_{\text{cell}} \to \mathbb{V}_{\text{cell}}^*\).

• The Preconditioner (\(P\)): The preconditioner is responsible for bridging the two spaces. It represents an operator that, given a dual vector \(\tilde{r} \in \mathbb{V}_{\text{cell}}^*\) (in cell-wise format), returns a primal vector \(z \in \mathbb{V}_{\text{cell}}\) that is an approximation of the solution satisfying the constraints defined by \(C\). \(P \colon \mathbb{V}_{\text{cell}}^* \to \mathbb{V}_{\text{cell}}\).

3.2 The Conjugate Gradient Algorithm↩︎

The Conjugate Gradient method naturally arises from the minimization of the quadratic energy functional \(\mathcal{J}(u) = \frac{1}{2} \langle A u, u \rangle - \langle \tilde{b}, u \rangle\) over the continuous primal space. The variation (gradient) of this functional yields the residual \(\tilde{r} = \tilde{b} - A u\). Because the operator \(A\) maps continuous fields to distributions, this residual resides in the dual space \(\mathbb{V}_{\text{cell}}^*\).

This poses a mathematical challenge: we cannot directly update the primal solution \(u\) using the dual residual \(\tilde{r}\), as adding a functional to a field is an ill-defined operation. The preconditioner \(P \colon \mathbb{V}_{\text{cell}}^* \to \mathbb{V}_{\text{cell}}\) serves as the mathematical bridge between these spaces. It maps the unassembled dual residual back to a continuous primal correction \(z = P(\tilde{r})\) while implicitly defining the preconditioned inner product structure which governs the Krylov subspace.

To minimize the energy functional along a continuous search direction \(p_k \in \mathbb{V}_{\text{cell}}\), the step size \(\alpha_k\) is obtained by enforcing orthogonality between the updated dual residual and the primal search direction (\(\langle \tilde{r}_{k+1}, p_k \rangle = 0\)). In standard PCG, this yields: \[\label{eq:alpha} \alpha_k = \frac{\langle \tilde{r}_k, z_k \rangle}{\langle \tilde{q}_k, p_k \rangle},\tag{6}\] where \(\tilde{q}_k = A p_k\). In the context of cell-wise storage, the numerator pairs a dual residual with a primal correction, and the denominator pairs a dual flux with a primal search direction. By Lemma 1, both of these duality pairings can be evaluated using only local, unassembled cell-wise operations.

In principle, the primal-dual structure above applies verbatim to standard PCG. Our motivation for going beyond it is hardware-driven: GPUs offer substantially higher throughput in reduced precision (single or even half), and the natural place to exploit this is inside the preconditioner [7], e.g.by running the multigrid V-cycle in lower precision or varying its cycle structure between iterations. Such a preconditioner is no longer a fixed linear operator. Standard PCG, however, relies on a static preconditioner to maintain the mutual \(A\)-conjugacy of the search directions; if \(P\) changes between iterations, its convergence degrades or fails outright. We therefore employ the Flexible Conjugate Gradient (FCG) method proposed by Notay [21]. Since the two methods differ only in the formula for \(\beta_k\), we state them in primal-dual form as a single Algorithm 1, with the two \(\beta_k\) variants given side by side: in both cases every scalar is a duality pairing of one dual and one primal vector, computable from local data by Lemma 1. Since the primal-dual reasoning is identical for PCG and FCG, we present the proofs for the more general FCG variant only.

FCG modifies the search direction update to explicitly enforce \(A\)-conjugacy against the varying preconditioner, removing the reliance on a static inner product geometry. The updated search direction is \(p_{k+1} = z_{k+1} + \beta_k p_k\), where: \[\label{eq:beta} \beta_k = \frac{- \langle \tilde{q}_k, z_{k+1} \rangle}{\langle \tilde{q}_k, p_k \rangle}.\tag{7}\] FCG preserves the necessary primal-dual alignment. The new preconditioned residual \(z_{k+1}\) is primal, and the local operator output \(\tilde{q}_k\) is dual. Thus, the FCG numerator \(\langle \tilde{q}_k, z_{k+1} \rangle\) computes the global \(A\)-conjugacy locally, circumventing the need to explicitly assemble the operator output.

This flexibility comes at the modest price of one additional scalar product per iteration. In standard PCG the numerator of \(\beta_k\), \(\langle \tilde{r}_{k+1}, z_{k+1} \rangle\), is reused as the numerator of \(\alpha_{k+1}\) in the next iteration, so \(\beta_k\) requires no new reduction. The FCG numerator \(\langle \tilde{q}_k, z_{k+1} \rangle\), by contrast, must be computed in addition to it, while the shared denominator \(\langle \tilde{q}_k, p_k \rangle\) is already available from the computation of \(\alpha_k\). Since each iteration is dominated by the operator and preconditioner applications, this extra local pairing is negligible in practice.

By adopting this primal-dual algorithmic structure, the entire Krylov solver is reduced to parallel, local cell-wise operations. The requirement for communication, interface accumulation, and hanging-node resolution is completely abstracted away from the Krylov subspace projection and delegated exclusively to the preconditioner. The resulting routine is summarized in Algorithm 1, covering both the static (PCG) and flexible (FCG) choices of \(\beta_k\). We establish that this local formulation is equivalent to solving the globally assembled system, by the theorem below.

As a direct corollary from the theorem above, provided the standard explicitly assembled FCG algorithm converges for the given problem, the communication-free algorithm minimizes the global energy functional and converges to the true continuous solution \(u = \mathcal{G} u^{\text{global}}\).

3.3 Simple Preconditioners and Direct Stiffness Summation (DSS)↩︎

As established in Algorithm 1, the communication-free FCG solver inherently requires a preconditioner \(P \colon \mathbb{V}_{\text{cell}}^* \to \mathbb{V}_{\text{cell}}\) to map the unassembled dual residual back to a continuous primal correction. A natural question arises: what constitutes the simplest possible, or trivial, preconditioner in this framework?

Algebraically, one might consider the identity matrix \(P = I\). However, directly equating a functional to a field without a change of basis is mathematically ill-posed. The identity operation between the dual space and the primal space is the Riesz isomorphism. Because the unassembled dual residuals in \(\mathbb{V}_{\text{cell}}^*\) represent partial, localized integrals, applying this mapping requires a topological transition.

To formalize this, we conceptually project the operations through the continuous global topology. First, applying the scatter operator \(\mathcal{G}^T \colon \mathbb{V}_{\text{cell}}^* \to \mathbb{V}_{CG}^*\) performs an additive assembly across shared element boundaries, reconstructing a complete, globally defined functional. Within this assembled space, the discrete Euclidean Riesz isomorphism \(\mathcal{I}_{\text{Riesz}} \colon \mathbb{V}_{CG}^* \to \mathbb{V}_{CG}\) acts algebraically as a numerical identity on the coefficient vectors, reinterpreting the accumulated integrals as continuous nodal values. Finally, the primal gather operator \(\mathcal{G} \colon \mathbb{V}_{CG} \to \mathbb{V}_{\text{cell}}\) extracts this assembled \(C^0\)-continuous field back into the redundant block format.

This composite mapping defines the operator \[\label{eq:dss95operator} S = \mathcal{G} \mathcal{I}_{\text{Riesz}} \mathcal{G}^T \colon \mathbb{V}_{\text{cell}}^* \to \mathbb{V}_{\text{cell}},\tag{9}\] which is known as Direct Stiffness Summation (DSS). The trivial preconditioner, which plays the role of the identity in the assembled setting, is given by \(P = S\). Note that \(S\) satisfies the continuity requirement of the equivalence theorem by construction: its image is exactly \(\mathcal{G}(\mathbb{V}_{CG}) = \text{Im}(C)\).

This composite formulation also extends to basic relaxation schemes, such as the Jacobi preconditioner, which requires an inverse diagonal modifier \(D_{CG}^{-1}\). In a classical assembled framework, scaling must occur in the global space before gathering. However, the diagonal scaling intertwines with the gather operator: let \(d \in \mathbb{V}_{CG}\) be the vector of diagonal entries of \(D_{CG}\), and let \(D_{\text{cell}} = \operatorname{diag}(\mathcal{G} d)\) be the diagonal operator on \(\mathbb{V}_{\text{cell}}\) whose entries are the gathered copies of \(d\). Since \(\mathcal{G}\) merely replicates entries, scaling before gathering equals gathering and then scaling each copy by the same factor, i.e.\(\mathcal{G} D_{CG}^{-1} = D_{\text{cell}}^{-1} \mathcal{G}\). Consequently, \[\mathcal{G} D_{CG}^{-1} \mathcal{I}_{\text{Riesz}} \mathcal{G}^T = D_{\text{cell}}^{-1} (\mathcal{G} \mathcal{I}_{\text{Riesz}} \mathcal{G}^T) = D_{\text{cell}}^{-1} S.\] This justifies evaluating the unscaled DSS first, and subsequently applying the pointwise scaling directly onto the localized cell vectors (i.e., \(P = D_{\text{cell}}^{-1} S\)).

The inter-cell communication is effectively moved into the preconditioner, while the outer Krylov solver remains local. While the operator \(S\) provides the essential bridge for Krylov preconditioners and the multigrid smoother, its definition explicitly relies on \(\mathcal{G}\) and \(\mathcal{G}^T\). Executing these sparse mappings translates to severe gather-scatter memory latency on modern hardware. The following section details how our algorithmic design allows for circumventing the global memory bottleneck.

4 Algorithmic Design of the DSS Operator↩︎

The objective of our algorithmic design is to eliminate the global gather-scatter operations. Rather than routing data through the indirect indices of a global vector, our DSS algorithms exploit the dense localization of the block layout: partial sums are exchanged directly across cell interfaces via symmetric point-to-point transmissions, reproducing the action of \(S\) while operating within the high-bandwidth, vectorized \(\mathbb{V}_{\text{cell}}\) storage format. The construction is layered in three stages of increasing topological generality, each presented in its own subsection: a dimensionally-split sweep that resolves all interfaces interior to a structured block (Section 4.1) and a set of kernels for the unstructured interfaces between macro-blocks (Section 4.2). The stages act on disjoint sets of degrees of freedom and compose into the full operator \(S\). A pre-/post-processing wrapper that further reduces non-conforming interfaces of \(h\)-adaptive meshes to the conforming case is deferred to Appendix 8.

4.1 Base Algorithm for Structured Grids↩︎

To achieve the high throughput required for modern GPU architectures, our approach decomposes the structured block interior summation into a sequence of dimensionally-split passes along the principal axes. In a 3D triangulation, a single degree of freedom at a cell corner is shared by \(2^3 = 8\) adjacent elements. While a traditional matrix-free implementation necessitates a complex gather-scatter operation to aggregate all eight contributions simultaneously, our multi-pass strategy achieves the same result through local one-to-one exchanges: a single pass along axis \(\alpha\) visits every interior interface normal to \(\alpha\), sums the two coincident copies of each interface DoF, and writes the sum back to both sides, so that after the pass both copies hold the assembled value of that face pair.

The key property is that vertex and edge DoFs are assembled as a byproduct of these face passes, without ever being addressed explicitly. Consider, in 2D, a vertex shared by four elements (Figure 2). Each element initially holds its own partial contribution. The \(y\)-pass sums the copies across the two horizontal interfaces, after which each copy holds the sum of its vertical pair of elements. The subsequent \(x\)-pass exchanges these already pairwise-summed values across the vertical interfaces, so each of the four copies ends up holding the total of all four contributions: the partial sums cascade around the corner. The same argument applies recursively in 3D: a vertex DoF shared by \(2^3 = 8\) elements receives its full sum after the \(z\), \(y\), and \(x\) passes, and an edge DoF shared by four elements after the two passes perpendicular to it. Sequentiality of the passes is essential. Each pass must observe the partial sums produced by the previous one, but only \(d\) kernel launches are required in total, independent of the polynomial degree and of the block size.

At runtime, we execute these passes using specialized kernels that permit fine-grained control over the memory hierarchy. The kernel maps thread blocks to batches of adjacent element interfaces. Because the block layout is uniform, extraction and injection are performed using statically unrolled loops and hard-coded face offsets, ensuring coalesced memory accesses. Furthermore, by processing one axis at a time and ensuring each interface pair is exclusively handled by exactly one thread, we safely update the boundary values without requiring hardware atomic operations. This translates a latency-bound graph-traversal problem into a throughput-bound streaming operation. By the end of this structured phase, all internal block degrees of freedom have reached their final assembled values, while the block boundary degrees of freedom contain partial sums representing their local block’s contribution.

4.2 Generalization to Block-Structured Grids↩︎

The dimensionally-split cascade of Section 4.1 relies on every interior vertex having the regular \(2^d\)-element neighborhood of a tensor-product grid, so that the axis passes enumerate exactly the right exchange partners. At the interfaces between macro-blocks this regularity is lost: a macro-vertex of an unstructured base grid may be shared by any number of blocks (e.g.three or five in 2D, Figure 4), and no sequence of axis-aligned pairwise exchanges visits all of its copies. These interfaces must therefore be resolved by dedicated kernels that traverse the macro-topology explicitly. Their cost is minor, as the affected DoFs form a lower-dimensional skeleton of the mesh, but their correct and divergence-free execution requires care.

To prevent race conditions and branch divergence on the GPU, the unstructured interfaces are partitioned by codimension into three mutually disjoint sets, each processed by its own kernel launch: Faces (macro-block quads, or in 2D the codimension-one segments, with their bounding lines and corners masked out), Lines (macro-block edges with their endpoints masked out, in 3D only), and Vertices (macro-block corners). The exclusive masking, where each kernel touches only the strict interior of its entity, iterating the open index range \(1,\dots,p-1\) along the entity and leaving the bounding sub-entities to the lower-dimensional passes, guarantees that the three launches write to non-overlapping DoFs, so they may run in any order without atomics or inter-launch synchronization. This inter-block phase operates on the macro-interface skeleton and is disjoint from the structured intra-block phase; it runs after the structured sweeps, consuming the partial block-boundary sums they produce and completing them into the globally assembled values.

Faces differ from the lower-dimensional entities in valence: a codimension-one interface is shared by exactly two blocks, so the Face kernel is a fixed pairwise exchange (cells \(a\) and \(b\)) rather than a variable-length gather. The only complication is orientation: unstructured mesh generators can connect the two blocks with arbitrary relative rotations or reflections, so before summing, side \(b\) must be re-indexed into the local frame of side \(a\). This is resolved on-the-fly with pure integer arithmetic, with no connectivity lookup: the relative orientation is packed into a small bitfield (in 3D three logical flags indicating an axis swap and flips along the first or second local coordinate) from which the matching index on side \(b\) is computed algebraically (Figure 4). The kernel then loads both copies, sums them, and stores the result back to each side under its own indexing.

Lines and Vertices, by contrast, carry a runtime valence: a macro-edge or macro-corner of an irregular base grid may be shared by any number of blocks. Their kernels gather the sharing set into a thread register, accumulate, and scatter the completed sum back to every copy; these are the only kernels in the entire DSS pipeline that consult a sharing list rather than a fixed-arity neighbor map. The sharing list is padded to a fixed width so the launch stays branch-free, and (like the Faces) each copy is re-indexed into a common canonical frame before accumulation. A one-dimensional reversal flag for lines, no orientation for vertices.

A subtlety distinguishes these entities from the structured and Face passes and is in fact the crux of the block-structured generalization. A macro-corner (or, in 3D, a macro-edge) is in general shared by several cells of the same block, and the structured pass has already summed those into a common per-block value that it replicates on every within-block copy. Adding all copies again across blocks would therefore multiply the result by the within-block multiplicity. The gather thus sums only one representative copy per block and scatters the total back to all copies, representatives and non-representatives alike. Concretely, take the midpoint \(X\) of a shared macro-edge in Figure 4: it is the corner at which two sub-cell faces meet, hence a Vertex entity, and in 2D it is shared by four cells (two in each of the two adjacent blocks). Initialize every copy to \(1\). The structured pass sums the two within-block copies on each side, leaving \(2\) on all four. A naive gather over all four copies would return \(2+2+2+2 = 8\); instead the Vertex kernel adds one representative per block, \(2 + 2 = 4\), and writes \(4\) to every copy, which is the correctly assembled value. The number of contributions to such an entity is therefore its number of sharing blocks, not its number of cell copies: two for the mid-edge node (a pairwise exchange), and three or five for the central macro-corner of Figure 4. The same reasoning applies verbatim to the 3D Line kernel, where a macro-edge interior is likewise reached by several cells of each block; the representative selection carries over unchanged.

5 Implementation↩︎

5.1 Memory layout: the BXYZE format↩︎

The field data reside in a multi-dimensional array with layout \([B,\, n_x,\, n_y,\, n_z,\, n_e]\) (which we refer to as the BXYZE layout), where \(B\) is a fixed block size, \((n_x, n_y, n_z)\) are the local DoF indices within an element, and \(n_e\) enumerates the elements of a block. Because the cell-wise storage decouples the layout from the global mesh topology, this ordering can be chosen freely to suit the hardware; the choice governs a trade-off between memory coalescing and element-contiguity, best seen at the two extremes. In an element-major layout (EXYZ), the element index is the fastest stride: threads operating on the same local DoF across consecutive elements read consecutive words: coalesced. However, operations that touch only a proper subset of an element’s DoFs (e.g.the DoFs of a single face during DSS or face restriction) find these fragments interleaved across distant elements, fragmenting the memory request pattern. Conversely, in a node-major layout (XYZE) all DoFs of one element are contiguous, so element-local sub-accesses are dense; but threads working on distinct elements stride by the full per-element footprint \((p+1)^d\), and coalescing across elements is lost.

The tension is thus between two classes of kernels. For volumetric operator evaluation, where every thread processes whole elements, EXYZ is the optimal layout: accesses are coalesced and nothing else matters. But the same layout is counterproductive for the DSS kernels of Section 4, which read and write face, edge, and vertex fragments: in EXYZ these fragments are scattered with element-sized strides, and the surface kernels’ throughput collapses. The blocked BXYZE layout is the compromise: grouping \(B\) spatially adjacent elements into one memory tile keeps accesses along the fastest axis coalesced for the volumetric kernels, while the per-element data remain compact enough that interface fragments span few cache lines. (Internal padding where fewer than \(B\) valid elements exist is a negligible concern: with \(B \le 128\), at most \(B\) wasted cells per block boundary are irrelevant at realistic problem sizes.) The layout also maps naturally onto sum factorisation: the 1D contractions operate independently along the \(B\) and \(n_e\) axes and translate directly into vectorised multi-dimensional array operations. We back these claims up empirically in Section 6.1, where the volumetric and DSS kernels are benchmarked across the layout design space.

5.2 Kernel implementation in Triton↩︎

Our entire solver stack is implemented in Triton [33], a domain-specific language and compiler for GPU kernels embedded in Python. For readers from the finite element community unfamiliar with it: Triton occupies a middle ground between high-level array frameworks (such as JAX/XLA [34]) and hand-written CUDA. A Triton kernel is written as a Python function operating on small multi-dimensional tiles (blocks of a tensor) rather than on individual threads; the programmer specifies what each tile-level program instance loads, computes, and stores, while the compiler decides the thread-level details: register allocation, vectorization of loads, shared-memory staging, and the mapping of tile operations onto the SIMT hardware. Notably, Triton bypasses the CUDA toolchain: kernels are compiled through its own MLIR-based pipeline directly to PTX (or, on AMD hardware, to the corresponding AMDGPU ISA), so no .cu source is ever generated, and the same kernel source runs on both NVIDIA and AMD GPUs.

This tile-based abstraction matches our cell-wise storage format exactly: a memory tile of the BXYZE layout is a natural Triton block, and the sum-factorized 1D contractions become tile-level tensor operations. The result is that kernels which would require hundreds of lines of carefully tuned CUDA are expressed in a few dozen lines of Python, while remaining competitive with vendor libraries in achieved bandwidth. Because Triton operates on PyTorch tensors, every kernel can be verified directly against a reference implementation written in high-level, abstract linear algebra on the same data, which is an invaluable property when developing the element kernels, where reproducing the sum-factorized tensor contractions by hand is highly error-prone, especially for the intricate tile decompositions of non-power-of-two polynomial degrees (\(p = 4, 5\)). We also note that a distributed, multi-GPU extension of the programming model exists [35], offering a natural path to scaling out the present single-device implementation.

Two implementation aspects deserve explicit mention. First, shared-memory bank conflicts: when a tile is staged in on-chip shared memory, which is divided into 32 banks, threads that simultaneously access different addresses in the same bank are serialized. Tensor-product index patterns are particularly prone to this, since the stride between consecutively accessed DoFs is often a multiple of the bank count. Triton’s compiler mitigates these conflicts automatically by swizzling (applying a bank-aware permutation to the shared-memory layout of staged tiles [33], [36]), so the contractions proceed without swizzling tricks customary in hand-written CUDA kernels (see for example [28]).

Second, power-of-two tile shapes: Triton requires every tile dimension to be a power of two. The per-element DoF count \((p+1)\) meets this only for \(p \in \{1, 3, 7\}\). For \(p = 2, 4, 5\), padding a single tile up to the next power of two would waste a substantial fraction of memory and bandwidth, so our data layout instead decomposes the DoF range of such elements into several tensors of power-of-two extent (e.g.\(5 = 4 + 1\)), and each kernel processes the resulting sub-tiles; this complicates the kernels but keeps the stored data dense. For \(p = 6\) we make the opposite choice and simply pad to \(8^3\): the utilization \((7/8)^3 \approx 0.67\) is still reasonable, and avoiding the multi-tensor decomposition removes its structural overhead. All supported degrees \(p = 1, \dots, 7\) are thus covered by one of the two strategies. This padding is the source of the gap between the theory and achieved TFLOP/s rates for \(p \ge 3\) (cf.Table 4): the hardware issues operations on the padded tiles, while the theoretical rate counts only the useful contractions.

Finally, the dominant arithmetic primitive (the 1D contraction of sum factorization) is expressed through Triton’s dot primitive, which the compiler lowers onto the GPU’s matrix units (tensor cores). The contractions of all elements in a tile share the same small 1D matrix (\(V\), or its transpose), so they batch naturally into the matrix-matrix shapes these units consume. The matrix units double the FP64 peak relative to the vector pipeline; it additionally frees the vector pipeline for address arithmetic and keeps the contraction off the critical path. At low polynomial degrees the 1D matrices are too small to fill the tensor-core tile shapes that Triton supported out of the box. I therefore implemented support for the small-tile MMAv2 instruction shape (\(8 \times 8 \times 4\)) in the Triton compiler, a contribution that has since been merged upstream.² On the benchmarked matrix-matrix multiplication shapes, Triton for some sizes slightly outperforms cuBLAS, which we attribute to the compiler’s software pipelining of the main loop, overlapping global-memory staging with tensor-core computation.

It is worth stressing, for readers from the finite element community, what the implementation does not contain. There is no global degree-of-freedom numbering and no element-to-global index map anywhere in the solver: every kernel reaches its operands through statically known offsets into the \([B,\, n_x,\, n_y,\, n_z,\, n_e]\) tiles. The only connectivity data that survive are the macro-block adjacency lists consulted by the unstructured DSS vertex- and line-kernels of Section 4.2, and these reference merely the lower-dimensional skeleton of the mesh. The gather-scatter index arrays that consume both storage and bandwidth in classical matrix-free frameworks thus have no counterpart in our data structures.

6 Numerical Results↩︎

Throughout this section we report throughput in giga-degrees-of-freedom per second (GDoF/s), counting each degree of freedom once. This normalizes the cost to the size of the problem actually being solved rather than to an implementation-dependent count of local element nodes, and it makes rates directly comparable across polynomial degrees and across operators. The unique-DoF count itself is obtained in one of two ways. On structured Cartesian grids it is known in closed form, \(N_\text{unique} = (n p + 1)^d\) for \(n\) cells per direction, and the reported rates are exact. On general (unstructured or \(h\)-adaptively refined) meshes we instead measure the cell-wise rate (total stored DoFs divided by apply time) and rescale it to unique DoFs by the large-mesh factor \(\big(p/(p+1)\big)^d\). This factor is the asymptotic interior ratio and ignores the surface DoFs of the domain boundary, which are shared by fewer elements; it therefore slightly overcounts the redundancy and undercounts the unique DoFs, so the resulting unique rates are a conservative lower bound on the true throughput.

One further caveat applies specifically to the DSS rates: for \(p > 1\) only the interface DoFs shared between elements take part in the direct stiffness summation (the interior DoFs of each element are never exchanged), so the per-unique-DoF DSS throughput is computed against the full DoF count even though a shrinking fraction of those DoFs is actually touched as \(p\) grows. This is deliberate: it keeps the DSS rate on the same per-DoF footing as the volumetric rates, so the two can be read on a common axis, and it is the reason the apparent DSS throughput rises with \(p\) even though the work per interface DoF does not.

6.1 Layout study: the volumetric–DSS trade-off↩︎

We begin by substantiating the central layout claim of Section 5.1: that the block size \(B\) of the \([B,\, n_x,\, n_y,\, n_z,\, n_e]\) layout governs an irreconcilable tension between volumetric operator evaluation and the surface DSS kernels, and that an intermediate \(B\) is the right compromise. The measurements are performed on a Cartesian grid obtained by uniform hierarchical refinement of a single coarse cell (\(16.7\) million active elements, \(\approx 1.07\) billion active DoFs at \(p=3\)); on this structured geometry every macro-interface is conforming and regular, so the reported DSS figures exercise the structured-path kernels of Section 4.1 alone, isolating them from the unstructured and hanging-node machinery. For each polynomial degree we sweep \(B\) over powers of two and, at every \(B\), retain the launch configuration that maximizes volumetric throughput, reporting the corresponding mass-operator volumetric and DSS rates per unique DoF.

Figure 6 (a) lays the trade-off bare. The solid curves (volumetric throughput) are essentially flat in \(B\): whole-element streaming is already coalesced once \(B \ge 2\), so enlarging the block buys nothing and, through padding and slightly worse cache behaviour, costs a few percent at the largest sizes. The dashed curves (DSS throughput) behave in the opposite way: at \(B=1\) each structured face exchange straddles a block boundary and the kernel is starved, but it climbs steeply and monotonically with \(B\) as more interfaces become intra-block and the axis sweeps run on dense, coalesced tiles, saturating only past \(B \approx 32\)–\(64\), an order-of-magnitude gain from the smallest to the largest block. The two curves cross around \(B \approx 16\)–\(32\): below it the surface kernels dominate the cost, above it they are essentially free while the volumetric rate has begun to erode. A block size of \(B = 16\)–\(32\) therefore captures almost all of the DSS speed-up at a negligible volumetric penalty, which is precisely the compromise adopted throughout the remainder of this work. Notably, the DSS rate improves with polynomial degree (the surface-to-volume ratio of the exchanged DoFs falls as \(p\) grows) so the trade-off only relaxes for the higher-order elements that are the main target of the method.

6.2 DSS throughput dependence on mesh structure↩︎

We now isolate the cost of the direct stiffness summation itself and ask how much the structured (dimension-split) path of Section 4.1 gains over the general unstructured path. Both are exercised on the same geometry (a unit cube under uniform refinement) so that the meshes are identical and only the code path differs: the structured path exploits the tensor-product regularity of the interfaces, whereas the unstructured path treats every shared DoF through the general indirection machinery. DSS carries essentially no arithmetic; each shared DoF is gathered, summed, and scattered. Its arithmetic intensity is therefore negligible and the kernel is memory bound, so we report only the achieved DRAM traffic and, alongside the GDoF/s rate, its speed-of-light (SoL): the fraction of peak memory bandwidth attained. Table 6 collects both metrics for \(p = 1, \dots, 7\).

Two trends stand out. First, the structured path sustains a high and roughly constant bandwidth utilization (around \(80\%\) of peak across all degrees), while the unstructured path plateaus markedly lower, near \(60\%\), the gap reflecting the scattered, less coalesced memory accesses of the general indirection. Second, the throughput advantage of the structured path widens with polynomial degree: at \(p = 1\) the two paths deliver essentially the same GDoF/s rate, but by the highest degrees the structured path is substantially faster. The \(p = 1\) tie is instructive. There the two paths run genuinely different algorithms: the structured pass streams every vertex of the Cartesian grid a fixed number of times by construction, which inflates its traffic and hence depresses its SoL even though it moves the same number of unique DoFs per second as the unstructured pass. Their equal throughput at \(p=1\) is thus a coincidence of two different cost structures meeting, not a sign that the paths are equivalent; as soon as \(p > 1\) and interior DoFs dilute the interface work, the structured path’s regularity pulls ahead.

Table 7 folds the two stages into the throughput of the entire assembled mass operator (one volumetric apply followed by a complete DSS, i.e.a single matrix–vector product of the continuous operator \(C^T A C\). The volumetric cost is fixed by the geometry (Cartesian or general/curved) and the DSS cost by the path (structured or unstructured), so the four rows bracket the realistic operating range. Throughput rises monotonically with \(p\), from \(3.0\,\)GDoF/s at \(p=1\), where the cell-wise redundancy is most punishing, to \(27.0\,\)GDoF/s at \(p=7\) on the Cartesian/structured path, as the surface-to-volume ratio of the exchanged DoFs falls and the coalesced volumetric apply comes to dominate the increasingly cheap DSS. The bottom row lists the libCEED BP1 mass-operator throughput [37] on the same A100 80 GB SXM device. Since libCEED’s BP1 is posed on a general (curved) geometry, the most directly comparable configuration is our general-geometry volume resolved through the structured DSS path; against it our entire operator is faster per unique DoF throughout, from \(1.6\times\) at \(p=1\) to \(4.4\times\) at \(p=5\), the margin widening with degree. This is the direct dividend of removing indirect gather–scatter from the operator: the volumetric apply streams coalesced cell-wise data and the DSS adds only the small, bandwidth-bound interface traffic of Table 6, whereas the assembled BP1 kernel pays the full uncoalesced gather–scatter cost on every apply.

6.2.1 DSS on unstructured meshes↩︎

The benchmarks so far ran on Cartesian or single-block geometries. We now exercise the DSS operator on meshes that are representative of production setups, reporting the DSS cost in isolation (the volumetric apply is excluded); the meshes are detailed in the captions and we read off the qualitative trends here.

The representative case is Figure 8: a curved ball domain obtained by \(6\) global refinements of an unstructured coarse grid, with no hanging nodes. Such a mesh is heavily block-structured in the interior of each coarse cell yet stitched together by genuine unstructured exchanges across the macro-block interfaces, which is the very mix a real application presents, and one that drives both the structured sweeps and the unstructured inter-block kernels at once. Its DSS throughput (Figure 8b) sits within a few percent of the structured-Cartesian ceiling of Table 6 across the whole degree range, well above the all-unstructured rate of the same table. Because the unstructured macro-interfaces form only a sparse skeleton of the mesh, they are slower per DoF, but touch too few DoFs to move the aggregate. Two facts reinforce each other here: a curved mapping leaves the assembly untouched (DSS moves only DoF values and references no geometric data, unlike the volumetric apply, whose cost the per-quadrature-point Jacobians inflate), so the general-mesh penalty of the combined operator lives in the volumetric stage; and a realistic, predominantly block-structured topology runs near the structured best case rather than the unstructured worst case.

7 Conclusions↩︎

In this work, the redundant, cell-wise vector is the persistent state of a high-order continuous finite element solver rather than a transient byproduct of matrix-free evaluation, and we have built the machinery that makes this practical on GPUs. The theoretical core is a primal-dual equivalence theorem: provided the preconditioner returns a continuous correction, the flexible conjugate gradient iteration executed on unassembled cell-wise data is identical, iterate by iterate, to the assembled solve, so inter-element communication is confined to a single direct stiffness summation inside the preconditioner and removed from the Krylov projection altogether. The algorithmic core is a DSS operator that realizes this summation without indirect addressing, atomics, or coloring: a dimensionally-split cascade of one-to-one face exchanges that accumulates edge and vertex sums as a byproduct of sequential axis passes, extended to unstructured macro-interfaces by disjoint topological kernels and to \(h\)-adaptive meshes by a shadow-cell wrapper that leaves the high-throughput sweeps untouched. Because cell-wise storage decouples the memory layout from the mesh topology, the layout itself becomes a tuning parameter, which we mapped across the blocked design space.

The payoff is measured on an A100. The volumetric operators run at the memory roofline: the mass operator transfers exactly the two doubles per DoF its formulation demands, and the Laplace operator, despite roughly three times the arithmetic, runs faster still, because its contractions fill the tensor cores more completely and leave it even more firmly memory-bound. The structured DSS sustains around \(80\%\) of peak bandwidth, and the entire assembled operator (volumetric apply plus a full DSS) exceeds the libCEED BP1 benchmark on the same hardware across the whole degree range, by up to roughly a factor of four, already at the lowest order. On a production-like curved mesh the assembly runs within a few percent of its structured ceiling, confirming that the unstructured machinery, though slower per DoF, is exercised on too thin a skeleton to matter. The complete stack is implemented in compact Triton kernels, runs unmodified on NVIDIA and AMD hardware, and contributed a small-tile tensor-core instruction shape upstream.

The price is memory footprint: the cell-wise format stores up to \(8\times\) as many values as the assembled vector at \(p=1\), a redundancy that falls steeply with degree and is mildest precisely in the high-order regime the method targets. The trade is a larger but coalesced footprint in exchange for eliminating the indirect gather-scatter traffic that bottlenecks classical pipelines, and the measurements show the exchange to be favourable at every degree tested.

Several directions follow naturally. The flexible Krylov variant was adopted precisely so that the preconditioner may change between iterations and, in particular, be evaluated in reduced precision [7]: because the solver is memory-bound, running the preconditioner in single or half precision directly and substantially lowers the bytes transferred per DoF (the quantity that governs performance throughout this work) at no cost to the accuracy of the outer iteration, which remains in double precision. The companion paper [32] develops the geometric multigrid that exploits this, and shows in addition that DSS across hanging nodes is not required by the solver, so the one bounded overhead identified here can be sidestepped. The framework is moreover agnostic to the smoother: overlapping patch smoothers of Schwarz type satisfy the consistency condition of the equivalence theorem by construction (Remark 2), each patch solve returning a continuous correction, so used as the preconditioner they render the outer iteration communication-free without any explicit DSS at all. The same agnosticism extends to the operator (any symmetric operator fits the theorem unchanged) opening the communication-free treatment to a far broader class of problems than the mass and Laplace operators benchmarked here.

7.0.0.1 Acknowledgements

The author would like to thank Mario Lezcano Casado (OpenAI) and Keren Zhou (George Mason University) for reviewing and accepting the author’s Triton patch enabling fp64 MMA support (\(8{\times}8{\times}4\) tensor-core instructions), on which the kernels of this work rely. Mario Lezcano Casado additionally fixed an independent Triton bug causing shared-memory bank conflicts for 64-bit operands, which substantially improved the throughput of the higher-degree kernels (\(p = 4, 5\)).

7.0.0.2 Declarations

Language models (Claude, Gemini) were used to produce unassembled initial drafts and to smooth the assembled text; as everywhere in this framework, assembly, and full accountability for all scientific content, remains with the author.

8 Extension of DSS algorithms onto Non-Conforming Interfaces↩︎

The construction below extends the DSS sweeps to \(h\)-adaptive meshes with hanging nodes. It currently works on structured (Cartesian) refinement; the unstructured hanging-node case is not yet robust, and the companion multigrid paper [32] shows that the solver does not require DSS across hanging nodes at all, so this path is not needed for the results of the present work. We record it here for completeness.

In the context of \(h\)-adaptive meshes with hanging nodes, interfaces may bridge elements at different hierarchical levels. To process these non-conforming boundaries without disrupting the highly optimized structured and unstructured kernels, we deploy a pre- and post-processing wrapper relying on Shadow Cells. A shadow cell is a structurally identical, conceptually dormant parent cell that is allocated in the cell-wise storage exactly like an active cell but exists exclusively to broker data across a refinement interface (Figure 9). Rather than treating hanging nodes as a special case during the summation sweeps, the algorithm restricts fine-level data to the coarse-level shadows, runs the standard structured and unstructured DSS on a topology that is globally conforming at the coarse level, and then prolongates the assembled data back to the fine cells. The non-conforming case is thereby reduced exactly to the conforming one, and the summation kernels remain completely unaware of adaptivity.

The pre-processing step executes a Restrict-Add and Zero-Out pattern. First, all shadow cells are explicitly zeroed to clear any potential stale data. Then, the dual residual degrees of freedom from the fine cells lying on the refinement edge are restricted to the coarse shadow cell using 1D tensor-product restriction operators \(\mathcal{R}_h\) and accumulated. Immediately following this accumulation, the corresponding values on the fine children are zeroed out. This zeroing process is critical to prevent double-counting when the subsequent standard DSS passes inevitably process the perpendicular axes.

Once the pre-processing is complete, the standard structured and unstructured DSS operations are executed, naturally including the populated shadow cells as valid participants in the global summation. In the post-processing phase, the assembled data in the shadow cells is lifted back to the fine grid. We inject the coarse shadow values back into the fine children using the prolongation operators \(\mathcal{P}_h\) and accumulate them. To prevent the injection from corrupting vertex values due to interpolation coupling, a save-restore pattern is used for shadow vertices before the face injection. Finally, the fine cells are zeroed out locally before receiving the accumulated values, and the shadow cells themselves are reset. This completely encapsulates the adaptivity logic outside of the high-throughput summation sweeps.

On a structured \(h\)-adaptively refined ball (\(5\) global plus \(4\) adaptive refinements, \(1.2\times10^{5}\) hanging interfaces) the wrapper’s throughput is reported in Table 8. The wrapper’s restriction and prolongation add work that grows as \((p{+}1)^2\), so the rate trails the conforming case—by a few tens of percent at low order, up to roughly a factor of two-and-a-half by \(p=5\)—yet it stays bounded and keeps rising with \(p\), since the non-conforming problem is reduced exactly to the conforming one and adaptivity never touches the high-throughput sweeps.

References↩︎