Matrix-Free Local Multigrid for $C^1$ Hermite Discretizations on GPUs

• The work extends a nested C1 multigrid setting toward cell-wise adaptive constraints and GPU execution, distinct from prior analyses of uniformly refined meshes.Source: §1 · p. 2
• The cell-wise primal-dual formulation enforces C1 conformity without forming the constraint operator, while geometric transfers supply both trace layers at refinement edges.Source: §1 · p. 3
• Even–odd operator evaluation and Chebyshev–Jacobi and vertex-block smoothers provide GPU-oriented alternatives within the same conforming multigrid hierarchy.Source: §1 · p. 3
• The numerical study is limited to the clamped biharmonic problem with cubic elements, although the construction is described as applicable to broader operator forms.Source: §1 · pp. 3–4
• Five or six iterations on uniform meshes, 85% of device peak memory bandwidth, and up to 4.23 million cells per second are reported in three-dimensional double precision on one NVIDIA A100.Source: §1 · p. 4

2 Biharmonic Problem and H2-conforming Discretization

• The model is the clamped biharmonic problem on a Cartesian domain, discretized in three dimensions on dyadically refined hexahedral meshes.Source: §2 · p. 4

2.1 Cell-Wise Representation and the Constraint Operator

• Cell-local degree-of-freedom copies represent the field, while the constraint operator characterizes admissible C1 vectors without requiring an assembled global vector.Source: §2.1 · pp. 4–5
• The constraint operator acts on the value and normal-derivative trace layers, enforcing matching-facet continuity and coarse-side interpolation at refinement edges without being formed.Source: §2.1 · p. 5

2.2 Tensor-Product Hermite Element

• Tensor-product Hermite basis layers encode the facet value and normal derivative, and geometric transfer scales reference derivative data by the child-to-parent width ratio.Source: §2.2 · pp. 5–6

3.1 Cell-Wise Inter-Grid Transfers

• Cell-wise prolongation evaluates the parent Hermite polynomial on each child, while restriction accumulates unassembled dual residual contributions from the children.Source: §3.1 · p. 6

3.2 Level-Wise Decomposition of an Adaptive Mesh

• Local multigrid decomposes the active-mesh residual by level and confines smoothing to cells covered at each level, keeping work proportional to active cells summed over levels.Source: §3.2 · pp. 6–7

3.3 Refinement-Edge Masking

• The smoother updates only interior degrees of freedom, while the volumetric operator reads refinement-edge data and the unmasked residual carries outward fluxes to coarser levels.Source: §3.3 · p. 7

3.4 The Local V-Cycle

• The local V-cycle accumulates restricted residuals and prolongates coarse corrections to set both value and normal-derivative edge data, with algebraic equivalence to the classical cycle after gathering.Source: §3.4 · pp. 7–8

4.1 Element Operator and Sum Factorization

• FDM regroups the six operator terms into two parts sharing first and last contractions, reducing evaluation from sum factorization’s 12 contractions to 10.The paper uses FDM as the baseline for measuring the decompositions.Source: §4.1 · p. 8, §4.2 · p. 9

4.2 Fast Diagonalization

• Fast diagonalization applies generalized-eigenvector transforms to two operator parts, whose correction terms reproduce the element operator exactly.The evaluation applies the operator, not its inverse.Source: §4.2 · pp. 8–9

4.3 Even–Odd Decompositions

• Signed reflection separates the cubic Hermite operator into parity sectors; splitting two axes yields four dense 16 × 16 products.The derivative degrees of freedom change sign under reflection, and the Cartesian operator commutes with signed reflection in each direction.Source: §4.3 · pp. 9–11

4.4 Comparison of Operator Evaluations

• For the cubic three-dimensional biharmonic operator, the two-axis decomposition had the highest measured throughput: 1.695 × 10^9 cells/s and 1735.4 GB/s, or 85% of peak memory bandwidth, 1.96× faster than FDM and 1.51× faster than dense matrix–matrix evaluation.Source: Table 2 · p. 12, §4.4 · p. 12

4.5 Roofline of the Operator Realizations

• The two-axis realization measures 3.5 FLOP/byte and 6.02 TFLOP/s, while FDM reaches 9.95 TFLOP/s at 11.3 FLOP/byte using tensor-core contractions.Roofline FLOP rates count issued arithmetic, including work the algorithm does not use.Source: §4.5 · pp. 12–13, Figure 1 · p. 12

5.1 Chebyshev–Jacobi Smoother

• Each Chebyshev–Jacobi stage applies a cell operator and pointwise update, then uses DSS to assemble the update and restore C1 continuity.DSS exchanges both Hermite trace layers across matching level interfaces; sequential axis passes propagate sums to edge and vertex copies without indirect addressing or atomics.Source: §5.1 · pp. 13–14, Figure 2 · p. 14
• Three Chebyshev stages smooth the upper spectral decade on [0.2, 2], leaving lower-frequency components to coarse-grid correction.The same construction uses 20 stages on the coarsest level instead of an exact solve.Source: §5.1 · p. 14

5.2 Vertex-Patch Smoother

• A vertex-patch update sums residual copies for the shared eight-DoF block, applies its dense inverse, and scatters the correction to every copy, preserving C1 conformity.On a uniform level, every interior vertex has the same cell configuration, so its patch matrix and inverse can be reused.Source: §5.2 · p. 15
• Eight parity-based colors allow concurrent non-overlapping patches within each color, while sequential colors make the sweep multiplicative and support symmetric pre- and post-smoothing.Forward pre-smoothing followed by reverse post-smoothing gives the symmetric composition needed for a conjugate-gradient preconditioner.Source: §5.2 · pp. 15–16

5.3 Performance of the Smoothers

• For n = 128, the patch sweep moves 8501 bytes per mesh cell, 8.30 times the field-data volume of one volumetric operator evaluation.The revisit factor is intrinsic to the colored multiplicative algorithm, separate from the hardware efficiency of its implementation.Source: §5.2 · p. 17, §5.3 · p. 17
• 1209 GB/s effective bandwidth is 59.3% of peak for the vertex-patch sweep, which takes 11.9 times as long as one volumetric operator evaluation.The bandwidth uses minimum field-data volume and may omit additional DRAM transactions for partial stores.Source: §5.3 · p. 18
• 1763 GB/s is reached by the addressing control, but its arithmetic and store volume both differ from the patch sweep.The control takes 14.281 ms versus 14.741 ms for the patch sweep, so this comparison does not isolate arithmetic or partial-store costs.Source: §5.3 · p. 18

6.1 Experimental Settings

• Experiments use cubic Hermite elements in three dimensions, FP64 arithmetic, and the device described in Section 4.4.The reported solve throughput is normalized by active cell count and reported alongside iteration count.Source: §6 · pp. 18–19, §6.1 · p. 19

6.2 Discretization and Constraint Verification

• The uniform-mesh L2 error rate is 3.97 to 4.00, and maximum jumps of both the value trace and physical normal derivative are zero.The reported measurements are for cubic Hermite elements on the unit cube.Source: §6.2 · p. 19, Table 5 · p. 20
• On uniform meshes, Chebyshev–Jacobi-preconditioned CG requires five or six iterations as the mesh grows from 512 to 2 097 152 cells.The symmetric patch cycle takes six or seven CG iterations; the stationary nonsymmetric patch cycle takes seven or eight.Source: §6.2 · p. 19, Table 5 · p. 20

6.3 Solve Throughput

• 4.23 Mcells/s is the finest-mesh throughput of the Chebyshev–Jacobi solver, versus 3.29 Mcells/s for the nonsymmetric patch cycle and 3.16 Mcells/s for its CG-accelerated symmetric counterpart.These are uniform-mesh solve throughputs.Source: §6.3 · p. 19, Table 6 · p. 20

6.4 Adaptive Refinement

• On locally refined meshes, errors remain close to and below those of the initial uniform meshes, with limited improvement under further refinement.The marked regions are fixed rather than selected by an error indicator; unrefined cells retain their mesh size and contribute to the error plateau.Source: §6.4 · pp. 21–22, Table 7 · p. 22
• The largest iteration-count increase occurs at the first local refinement; across tested adaptive meshes, the two preconditioners’ counts differ by at most one.For nuni = 5, Chebyshev–Jacobi counts after the first refinement are seven for the octant and eleven for the ball, while the ring count falls from sixteen to thirteen at later depths.Source: §6.4 · p. 22, Table 7 · p. 22

6.5 Refinement-Edge Geometry and Iteration Counts

• Iteration counts rise with the fraction of level cells incident to refinement edges: at the first adaptive step, 9%, 20%, and 74% edge-cell shares correspond to increases of one, five, and ten Chebyshev–Jacobi iterations.The ordering persists across both initial resolutions and all tested adaptive depths.Source: §6.5 · p. 23, Table 8 · p. 23

7 Conclusion

• The cell-wise Hermite formulation yields the expected fourth-order L2 convergence, with vanishing measured jumps in value and normal derivative across matching interfaces.Source: §7 · p. 23
• At refinement edges, the formulation constrains both trace layers and scales reference derivatives under geometric prolongation.Source: §7 · p. 23
• 85% of peak memory bandwidth is reached by the two-axis even–odd decomposition, which splits the element matrix into four dense 16 × 16 sector blocks.It has the highest throughput among the tested operator evaluations.Source: §7 · p. 23
• On the finest uniform mesh, Chebyshev–Jacobi preconditioned conjugate gradients reach 4.23 Mcells/s in five iterations, 18 times Witte et al.’s double-precision rate at the same degree.The symmetric patch preconditioner takes six iterations and reaches 3.16 Mcells/s; its sweep requests 8.3 times the field-data volume of one volumetric operator evaluation.Source: §7 · p. 23
• On adaptive meshes, the two preconditioners’ iteration counts differ by at most one; counts vary with refinement region, but the interface-area relationship is not isolated from other geometric differences.On the largest octant mesh, Chebyshev–Jacobi and patch preconditioning reach 3.56 and 2.77 Mcells/s, respectively.Source: §7 · p. 24

Use of generative AI

• The author used Claude (Anthropic) to assist with drafting and editing, reviewed and edited all content, and takes full responsibility for the publication.Source: Use of generative AI · p. 24