A Novel Matrix-Free Approximate Elemental Inverse for Efficient Preconditioning of Discontinuous Galerkin Spectral Elements

• Matrix-free evaluation avoids storing assembled elemental or global sparse matrices, while sum-factorisation reduces work through a sequence of one-dimensional tensor contractions.These approaches reduce memory traffic and floating-point operations, matching the high arithmetic intensity suited to modern hardware.Source: §1 · p. 1
• Efficient preconditioning remains a challenge for DG spectral-element Poisson solvers, which lack the mass term that makes implicit viscous Helmholtz systems more diagonally dominant.Hybrid multigrid and low-order-refined auxiliary-space preconditioning are identified as promising approaches, though the latter still requires an assembled auxiliary system.Source: §1 · p. 2
• Point Jacobi can be insufficient for strongly coupled modal bases, whereas dense block-Jacobi inverses improve local coupling treatment at substantial high-order storage and application cost.Earlier work also explored inexact local solves, transformed patch operators, and tensor-product approximations to reduce local-solve costs.Source: §1 · p. 2
• Extending tensor-product local solvers to collapsed-coordinate modal elements is difficult because their basis functions have coupled index ranges.The paper seeks a lower-cost local inverse that can compensate for additional solver iterations.Source: §1 · p. 3
• The proposed approximate elemental inverse combines inverse modal-to-nodal transforms with an element-dependent sparse matrix and uses sum-factorisation for efficient application across several element shapes.It is evaluated both as a standalone preconditioner and as a multigrid smoother; the largest gains are reported for high-order hexahedra and prisms.Source: §1 · p. 3

• The model is the Poisson equation with prescribed Dirichlet and Neumann boundary data, discretised using shape-dependent modal expansions on mapped elements.The discontinuous spectral-element space uses order-p expansions; element interfaces may carry different data on their two sides.Source: §2 · pp. 3–4
• SIPG imposes interface coupling through symmetric consistency, average-gradient flux, and jump-penalty surface terms, with boundary conditions imposed weakly.The penalty depends on order and element geometry through a local trace estimate; the largest adjacent-element value is selected for symmetric stability.Source: §2 · p. 4
• The discrete SIPG matrix is symmetric positive definite under the stated assumptions, enabling solution by conjugate gradients with convergence dependent on system conditioning and preconditioner effectiveness.Source: §2 · p. 4
• The SIPG Poisson operator is implemented in Nektar++ using matrix-free evaluation and sum-factorisation to reduce memory footprint and improve cache locality and performance.CPU SIMD instructions are used with a cache-friendly data layout.Source: §2 · p. 4

3 Polynomial multigrid preconditioning

• The framework uses p-multigrid for unstructured high-order discretisations because it avoids mesh agglomeration, with the standard preconditioned conjugate-gradient method requiring a fixed, linear, SPD preconditioner.Source: §3 · p. 4

3.1 Coarse-level solver

• Polynomial multigrid uses p=1 as its coarsest level, with Chebyshev iteration preconditioned by AMG to define a fixed linear SPD coarse solve.The low-order coarse system is assembled, with its cost considered acceptable because it is much smaller than the original high-order problem.Source: §3.1 · p. 5

3.2 Restriction and prolongation

• Restriction is the transpose of prolongation, which maps coefficients from a coarse polynomial level to the finer level.Prolongation can use L2 projection or basis interpolation; for nested spaces with full-column-rank basis matrices, interpolation gives the same result as L2 projection and shares matrices by element type.Source: §3.2 · p. 5

3.3 Smoother choices

• Chebyshev iterations serve as pre- and post-smoothers, while their effectiveness and cost depend strongly on the local preconditioner.At high order, smoothing-step cost matters because convergence analyses require the number of steps to grow with polynomial degree; point Jacobi may be insufficient for strongly coupled modal bases, while dense block Jacobi is costly.Source: §3.3 · p. 5

4.1 Modal-to-nodal transform and sparse inverse approximation

• The preconditioner uses a reference-basis left transform shared across like elements and a symmetric sparse matrix carrying element-specific geometry and local operator information.The transform is motivated by the clearer near-diagonal structure of element blocks in nodal bases and the spatial meaning of nodal degrees of freedom.Source: §4.1 · p. 6

4.1.1 Projection space

• Projection spaces include tensorial GLL grids, the original quadrature grid, and non-tensorial nodal point sets; collapsed-coordinate grids can contain more points than modal coefficients.The construction seeks a full-column-rank basis matrix and a left inverse; when the grid is overdetermined, the Moore–Penrose pseudoinverse is one option.Source: §4.1.1 · p. 6

4.1.2 Sparse patterns

• The tested sparse patterns range from diagonal Band<N> matrices to nearest-point and first-layer stencils, trading storage and computation against retained couplings.The first-layer stencil is largest, while practical patterns must remain sparse enough to offset transform costs.Source: §4.1.2 · pp. 6–7

4.2 Sum-factorisation for inverse transforms

• Applying a dense left inverse with GEMM would make transform work more costly than the sparse matrix application, motivating sum-factorised inverse transforms.Source: §4.2 · p. 7

4.2.1 Tensorial and partial tensorial sum-factorisation

• Hexahedral transforms use tensor-product contractions, while prismatic transforms partially factor the triangular basis and extrusion direction; inverse or pseudoinverse factors can be applied correspondingly.The paper names the hexahedral constructions SumFacInverse and SumFacPseudoInv.Source: §4.2.1 · p. 7

4.2.2 Generalised sum-factorisation in collapsed coordinates

• For tetrahedra and pyramids, generalised sum-factorisation handles collapsed-coordinate bases by reversing conditional contractions and correcting the modes absent from the uncorrected basis.The resulting SumFacLeftInv is an exact left inverse when the correction matrix is nonsingular, and the same principle extends to prisms.Source: §4.2.2 · p. 8

4.3 Implementation and performance optimisation

• The element-dependent sparse matrix stores O(P^3) entries per element and can have low arithmetic intensity, motivating interleaved SIMD and pattern-specific kernels.Element vectors and temporary workspace also require O(P^3) storage.Source: §4.3 · p. 9
• Sum-factorisation makes the approximate inverse more cost-effective than the exact elemental inverse by reducing inverse-transform work below the O(P^6) dense cost.Inverse-transform costs are O(P^4) for hexahedra, tetrahedra and pyramids, and O(P^5) for prisms; the current setup still constructs dense elemental inverses to obtain sparse entries.Source: §4.3 · pp. 8–9

4.4.1 Energy-scaled least-squares problem

• LSFit energy-scales the inverse-fit error so the exact inverse maps to identity, then solves for symmetric sparse entries with diagonally preconditioned CGLS.The iterative solve starts from the truncation-based sparse matrix and uses at most 30 iterations.Source: §4.4.1 · pp. 9–10
• The energy-scaled objective penalizes eigenvalue deviations from unity, but does not directly minimize condition number or guarantee faster global convergence.A smaller fitting residual alone therefore does not establish fewer solver iterations.Source: §4.4.1 · p. 9

4.4.2 Approximation quality and iterative convergence

• The normalized residual is the RMS deviation of local preconditioned eigenvalues from unity; a global condition-number bound requires every block spectrum to lie in a common positive interval.The least-squares objective measures combined local eigenvalue error, not the bound itself, so its convergence effect must be assessed in the complete solver and real cases.Source: §4.4.2 · p. 10
• Symmetry does not ensure positive definiteness, so diagonal entries of the sparse matrix may need correction, which can increase the fitting residual.The paper describes an empirical correction procedure in Appendix A.3.Source: §4.4.2 · p. 10

5.1 Benchmark setup

• Benchmarks solve SIPG Poisson problems on regular and twisted structured meshes; hexahedral cells are split into 2 prisms, 6 tetrahedra, or 6 pyramids, with normalized-residual tolerance fixed at 10^-9.Reported solver timings exclude setup time.Source: §5.1 · pp. 10–12

5.2 Impact of penalty factor and basis type on convergence

• For the tested shapes, the optimal penalty scale is generally about 0.2–0.25, except p = 1 tetrahedra at about 0.32; the tests use 0.32 on regular and 0.4 on distorted meshes.With point Jacobi, nodal hexahedral bases converge much faster than modal bases, especially at high orders; block-Jacobi iteration counts are similar for both bases.Source: §5.2 · p. 12

5.3 Selecting the fastest approximate elemental inverse preconditioner

• Tensorial GLL projection gives the lowest iteration counts among the tested projection spaces, while inverse-transform choices depend on element shape: SumFacPseudoInv is preferred for prisms and SumFacLeftInv is the available sum-factorised option for tetrahedra and pyramids.For tetrahedra and pyramids, SumFacLeftInv requires more iterations than FullPseudoInv but still yields shorter solver time; in hexahedra, SumFacPseudoInv is equivalent to FullPseudoInv at lower cost.Source: §5.3 · p. 13, Table 3 · p. 13, Table 4 · p. 14
• LeastSquaresFit particularly helps recover convergence for medium-to-high-order tetrahedra and pyramids with SumFacLeftInv, but can raise setup time by more than tenfold in tetrahedral cases.For hexahedra and prisms using SumFacPseudoInv, fitting does not consistently improve convergence; setup overhead may be amortised over repeated solves or reduced by reusing results for identical local operators.Source: §5.3 · pp. 13–14, Table 5 · p. 14, Table 6 · p. 14
• Sparse-pattern effects depend on shape and order: Band1 is fastest for the shown hexahedral cases, while NearestPnt is fastest for the shown high-order prism, tetrahedral, and pyramidal cases.For the shown medium- or low-order prisms and p = 2 tetrahedral and pyramidal cases, Band1 is fastest; denser patterns can reduce iterations but may cost more to apply.Source: §5.3 · pp. 14–15, Table 7 · p. 15

5.4 Applying the approximate elemental inverse to multigrid smoothing

• On regular meshes, the approximate inverse smoother has iteration counts similar to block Jacobi but is cheaper to apply; its time advantage is largest for high-order hexahedra and prisms, reaching an order of magnitude for p = 8 hexahedra.For low-order cases, the advantage over block Jacobi is modest.Source: §5.4 · p. 16, Table 8 · p. 16
• On twisted meshes, the approximate inverse remains faster than block Jacobi for hexahedra and prisms, but its advantage shrinks; results for tetrahedra and pyramids are more distortion-sensitive, including a p = 8 tetrahedral outlier with substantially more iterations than block Jacobi.The authors caution that it is not guaranteed to be the best choice for collapsed-coordinate elements on distorted meshes.Source: §5.4 · p. 16, Table 8 · p. 16

5.5 Comparison with the continuous Galerkin solver

• SIPG with p-MG and the approximate inverse is faster than continuous Galerkin on medium-to-high-order regular hexahedral, prismatic, and pyramidal meshes, while continuous Galerkin is faster on tetrahedral and low-order workloads.With a low-energy block preconditioner, continuous Galerkin is faster than SIPG in most twisted-mesh cases.Source: §5.5 · p. 17, Table 9 · p. 17

6 Conclusions

• The inverse combines shared inverse modal-to-nodal transforms with an element-dependent sparse matrix, using sum-factorisation variants across hexahedral, prismatic, tetrahedral and pyramidal elements.The structure reduces storage and memory access compared with the exact element-block inverse.Source: §6 · p. 17
• The approximate inverse generally needs more iterations than the exact block inverse, but its lower application cost can shorten solver time, with the greatest gains on hexahedral and prismatic meshes in standalone and multigrid use.The gains on tetrahedral and pyramidal meshes are smaller and more configuration-dependent.Source: §6 · pp. 17–18
• On tetrahedral and pyramidal meshes, least-squares fitting helps recover approximation quality, but lower fitting residuals do not ensure faster global convergence or positive definiteness.The approximate inverse is not guaranteed to outperform the exact elemental inverse across all workloads; setup cost can be amortised over repeated solves.Source: §6 · p. 18
• On medium-to-high-order regular hexahedral, prismatic and pyramidal meshes, SIPG can match or beat continuous Galerkin solver times, but continuous Galerkin remains faster for tetrahedral and low-order workloads and most distorted meshes.SIPG is particularly harder to solve on twisted tetrahedral and pyramidal meshes.Source: §6 · p. 18
• The matrix-free DG smoother remains useful in Fehn et al.'s hybrid multigrid framework, which uses continuous-space coarse-grid correction while retaining the discontinuous problem at the finest level [17].Finest-level smoothing usually accounts for most of the computational cost.Source: §6 · p. 18

A.1 Diagonal preconditioning and normal equations

• The diagonal entries of the least-squares normal matrix are squared Frobenius norms of the corresponding column outer products, and the unknowns are scaled as x_ij = √m_ij V_ij.Only pairs with m_ij > 0 are active; zero-scale pairs are omitted.Source: §A.1 · p. 18
• With the scaled unknowns, the symmetric sparse-pattern matrix V(x) defines F(x) = TV(x)T^T, and the least-squares objective is 1/2 ∥I − F(x)∥_F^2.The adjoint and forward operators use matrix products and restriction to the sparse pattern, avoiding assembly of the design matrix or normal matrix.Source: §A.1 · p. 18

A.2 Conjugate-gradient iteration

• The algorithm takes projected entries, active scales, the operators, a relative tolerance and an iteration limit as inputs, and returns an approximate least-squares solution.Source: §A.2 · p. 19
• Each iteration applies F to the search direction and F* to the updated matrix residual, then updates the solution and search direction by conjugate-gradient recurrences.The algorithm stops when the residual criterion is met or the iteration limit is reached.Source: §A.2 · p. 19
• CGLS solves the scaled normal equations using only F and F*, with matrix residual R_k = I − F(x_k) and normal-equation residual g_k = F*(R_k).The diagonal preconditioning is incorporated through the scaled unknowns.Source: §A.2 · p. 19

A.3 Positive-definiteness correction

• After fitting, the eigenvalues of the transformed matrix are checked because the least-squares iteration preserves symmetry but does not ensure positive definiteness.Source: §A.3 · p. 19
• If the eigenvalue condition fails, a positive shift is doubled from 10^-6 until the condition is satisfied; this empirical correction can increase fitting residual and does not exactly minimise the least-squares objective.Source: §A.3 · p. 19

• Figure 5 compares conjugate-gradient iterations and L2 error across penalty scale factors for hexahedral, prismatic, tetrahedral, and pyramidal meshes.Source: Figure 5 · p. 20
• Together, the figures let readers compare how penalty-scale variation is represented for regular versus twisted meshes using both convergence and accuracy measures.Source: Figure 5 · p. 20, Figure 6 · p. 20
• Figure 6 presents the same iteration and L2-error comparison across penalty scale factors for twisted hexahedral, prismatic, tetrahedral, and pyramidal meshes.Source: Figure 6 · p. 20

Regular-mesh preconditioner configurations

• For regular hexahedral, prismatic, tetrahedral, and pyramidal meshes, Figures 7–10 report iteration counts and solver times across preconditioner configurations.The approximate inverse uses a tensorial GLL projection space in Figures 7, 8, 9, and 10.Source: Figure 7 · p. 21, Figure 8 · p. 22, Figure 9 · p. 23

SIPG versus continuous Galerkin

• Figures 11 and 12 compare SIPG and continuous Galerkin solver iterations and times on regular and twisted hexahedral and prismatic meshes.The SIPG solver uses p-MG, with labels identifying preconditioners inside the Chebyshev smoother.Source: Figure 11 · p. 25, Figure 12 · p. 26
• Figures 13 and 14 make the same SIPG–continuous Galerkin comparison for regular and twisted tetrahedral and pyramidal meshes.The SIPG solver uses p-MG, with labels identifying preconditioners inside the Chebyshev smoother.Source: Figure 13 · p. 27, Figure 14 · p. 28