remarkRemark hypothesisHypothesis claimClaim A matrix-free ILU
realizationDaniel Drzisga, Andreas Wagner, Barbara Wohlmuth
scaletikzpicturetowidth[1] A matrix-free ILU realization based on
surrogatesDaniel Drzisga, Andreas Wagner, Barbara Wohlmuth

# A matrix-free ILU realization based on surrogates

Daniel Drzisga22footnotemark: 2 Andreas Wagner22footnotemark: 2 Corresponding
author. This work was partly supported by the German Research Foundation
through grant WO671/11-1. Barbara Wohlmuth Lehrstuhl für Numerische
Mathematik, Fakultät für Mathematik (M2), Technische Universität München,
Garching bei München<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Matrix-free techniques play an increasingly important role in large-scale
simulations. Schur complement techniques and massively parallel multigrid
solvers for second-order elliptic partial differential equations can
significantly benefit from reduced memory traffic and consumption. The matrix-
free approach often restricts solver components to purely local operations,
for instance, to the most basic schemes like Jacobi- or Gauss–Seidel-Smoothers
in multigrid methods. An incomplete LU (ILU) decomposition cannot be
calculated from local information and is therefore not amenable to an on-the-
fly computation which is typically needed for matrix-free calculations. It
generally requires the storage and factorization of a sparse matrix which
contradicts the low memory requirements in large scale scenarios. In this
work, we propose a matrix-free ILU realization. More precisely, we introduce a
memory-efficient, matrix-free ILU(0)-Smoother component for low-order
conforming finite-elements on tetrahedral hybrid grids. Hybrid-grids consist
of an unstructured macro-mesh which is subdivided into a structured micro-
mesh. The ILU(0) is used for degrees-of-freedom assigned to the interior of
macro-tetrahedra. This ILU(0)-Smoother can be used for the efficient matrix-
free evaluation of the Steklov–Poincaré operator from domain-decomposition
methods. After introducing and formally defining our smoother, we investigate
its performance on refined macro-tetrahedra. Secondly, the ILU(0)-Smoother on
the macro-tetrahedrons is implemented via surrogate matrix polynomials in
conjunction with a fast on-the-fly evaluation scheme resulting in an efficient
matrix-free algorithm. The polynomial coefficients are obtained by solving a
least-squares problem on a small part of the factorized ILU(0) matrices to
stay memory efficient. The convergence rates of this smoother with respect to
the polynomial order are thoroughly studied.

###### keywords:

ILU-Smoother, multigrid, hybrid grids, polynomial surrogates, matrix-free

65F55, 65N55

##  1 Introduction

The incomplete LU(0)-factorization[34] (ILU) approximates an LU factorization
by retaining the sparsity pattern of the original matrix. For strongly
anisotropic problems in 2D, it is often used as a smoother within multigrid
algorithms since its convergence rates are more stable than the ones of
simpler smoothers like the Gauss–Seidel- or Jacobi-Smoothers[47, Sec. 7.8].
This property carries on to anisotropic 3D problems in which the coupling in
one spatial direction is dominant while other schemes have to be used if two
of the spatial directions are dominant[26]. The related thresholded ILU-
Smoother was recently used for p-multigrid in isogeometric analysis[41, 42] or
as a smoother for the wave equation[45].

Besides its usage as a smoother, incomplete factorizations like the ILU are
used as preconditioners[5, 40], for instance in problems involving the
incompressible Stokes equation[24] or in electromagnetic scattering[31]. An
algorithm for a communication avoiding ILU(0) preconditioner in the high-
performance context was introduced in [22]. Algorithms for the efficient
parallel assembly of thresholded ILU preconditioners can be found in [3]
including adaptions to GPUs in [4, 30].

Matrix-free methods are becoming increasingly prevalent within finite-element
frameworks[29, 46, 27]. For instance, large scale mantle-convection
simulations typically operate on scales on which storing the discretization
matrices is not always feasible [6]. On the other hand, reducing the memory
traffic by not requiring to load a matrix from memory has the potential to
result in faster algorithms on today’s hardware. This generates interest in
adapting old matrix-based algorithms to the matrix-free context. For non-local
factorization algorithms like the ILU, this poses a tremendous challenge as
the matrix entries cannot be locally computed on-the-fly.

For structured grids, several techniques exist to approximate matrices for an
efficient evaluation. For instance, stencil-scaling techniques that work for
both scalar[7] and vectorial[19] equations. However, since the ILU-approach
relies on a matrix factorization which cannot be computed locally, these
approaches are inapplicable. In our work, we propose a matrix-free ILU
realization on structured subgrids based on surrogates. Here, the discrete
matrix, which usually approximates a continuous operator is additionally
approximated by surrogate polynomials [10, 9, 17, 8, 18]. These techniques can
also be adapted for hybrid structured grids which are extensively used in [11,
12, 13, 27] and consist of a coarse unstructured macro-grid which is
subdivided into a fine structured micro-grid. The former gives the approach
enough flexibility to represent relevant domains while the latter provides the
computational advantages of structured grids.

In this work, we apply the surrogate methodology to our factorized ILU matrix.
In the interior of the highly structured grids, we utilize an ILU
factorization and approximate the resulting matrix by surrogate polynomials.
This approximation is formed in a memory-efficient way such that the memory
costs stay within sensible bounds. We therefore obtain an efficient solver in
the interior of our structured grid.

To illustrate the potential of our approach, we provide two examples of how
the matrix-free ILU can be used on hybrid grids: Our main application is the
approximation of the Steklov–Poincaré operator for the Laplacian in a matrix-
free way. This operator is a main ingredient of many non-overlapping domain-
decomposition methods and therefore efficient algorithms for its evaluation
are highly relevant, see [15, 28, 43, 38, 32] and references therein. It
formally requires the exact inversion of an elliptic equation inside a
subdomain for which a multigrid method can be efficiently applied. By using
the ILU-factorization as a smoother within this inner multigrid, the inversion
becomes robust with respect to distortions along one axis. In the
supplementary material a second application is provided in which we extend the
subgrid ILU-Smoother to a smoother on the global grid.

The article is structured as follows: In Section 2, we describe the problem,
introduce the notation and present the Steklov–Poincaré operator. In Section
3, we introduce an ILU formulation that is amenable to a matrix-free
algorithm. Next, we introduce a reordering strategy on our hybrid mesh, to
optimize its performance as a smoother inside single subdomains. Finally, we
introduce the matrix-free surrogate ILU in Section 4 and compare its
asymptotic convergence rates within a multigrid algorithm to the matrix-based
ILU. We conclude with a short outlook and summary in Section 5.

##  2 Hybrid grids

In this section, we describe our model problem in the context of a low-order
conforming finite element discretization on hybrid grids. Hybrid grids combine
the flexibility of unstructured grids with the computational advantages of
structured grids [11, 12, 13, 33]. In addition, they provide a natural domain
partitioning that can be used to distribute the work to different nodes.

###  2.1 Preliminaries and notation

In the weak form of a Poisson-type equation, $-\operatorname{div}(K\nabla
u)=f$, on an open domain $\Omega\operatorname{\subseteq}\mathbb{R}^{3}$ with
homogeneous Dirichlet boundary conditions on
$\Gamma_{D}\operatorname{\subseteq}\partial\Omega$, natural boundary
conditions on $\partial\Omega\setminus\Gamma_{D}$ and an inhomogeneous,
bounded, symmetric, uniformly positive-definite diffusion tensor
$K(x):\Omega\to\mathbb{R}^{3\times 3}$, we obtain the bilinear form

(1) $a(u,v)=\int_{\Omega}\nabla u(x)^{\top}K(x)\hskip 0.25pt\nabla
v(x)\,\mathrm{d}x,\quad u,v\in V=\left\\{u\in
H^{1}(\Omega)\,:\,u|_{\Gamma_{D}}=0\right\\}.$

This includes the special case of a bounded, uniformly-positive scalar
material parameter $\kappa:\Omega\to\mathbb{R}$ by setting
$K=\kappa\,\textmd{Id}_{3}$, where $\textmd{Id}_{3}\in\mathbb{R}^{3\times 3}$
is the identity matrix. One application of the full diffusion tensor, would be
the pull-back of a blending function which maps a simple tetrahedral domain to
a more complex domain, thereby providing a better approximation of the domain
boundary. Given a load $f\in L^{2}(\Omega)$ which defines the linear form
$F(v)=\int_{\Omega}f\hskip 0.25ptv\,\mathrm{d}x$, we obtain the standard
variational problem: $\text{Find }u\in V\text{ satisfying }a(u,v)=F(v)\text{
for }v\in V.$

The typical approach in HHG[11, 12, 13] and HyTeG[27] is to discretize the
domain $\Omega$ with a coarse, possibly unstructured, simplicial
triangulation. This so-called macro-mesh consists of macro-vertices
$\mathcal{V}_{H}$, macro-edges $\mathcal{E}_{H}$, macro-faces
$\mathcal{F}_{H}$ and macro-tetrahedra $\mathcal{T}_{H}$. All macro-primitives
are referred to as
$\mathcal{P}_{H}=\mathcal{V}_{H}\cup\mathcal{E}_{H}\cup\mathcal{F}_{H}\cup\mathcal{T}_{H}$.
Based on this initial grid, we construct a hierarchy of $L\in\mathbb{N}$,
grids $\mathcal{T}=\\{\mathcal{T}_{h_{l}},\;h_{l}=2^{-l}H,\;l=2,\dots L+1\\}$
by successive global uniform refinement. The choice to start in the multigrid
hierarchy with $l=2$ guarantees that each macro-element contains at least one
interior element, which simplifies the notation in our algorithms. As it is
standard, each of these refinements is achieved by subdividing all elements in
3D into $8$ sub-elements. For details of the refinement in 3D, we refer to
[14]. Due to this refinement process, the element neighborhood at each vertex
in the interior of a macro element is always the same. The whole process of
the hybrid grid mesh setup is schematically depicted in Figure 1.

Associated with $\mathcal{T}_{h_{l}}$, is the space
$V_{h_{l}}\operatorname{\subseteq}V$ of piecewise linear conforming finite
elements

$V_{h_{l}}=\\{v\in V:v|_{t}\in\mathcal{P}_{1}(t)\text{ for each
}t\in\mathcal{T}_{h_{l}}\\}.$

Figure 1:  Hybrid-grid refinement procedure in 3D for a clipped tetrahedron in
a cubic macro-mesh. The DoF belonging to the index sets $\mathcal{I}_{v}$,
$\mathcal{I}_{e}$, $\mathcal{I}_{f}$, $\mathcal{I}_{t}$ and
$\mathcal{I}_{\partial t}$ are illustrated by different colors and shapes.

Let $\phi_{i}\in V_{h_{l}}$ and $\phi_{j}\in V_{h_{l}}$ be the scalar-valued
linear nodal basis functions associated with the $i$-th and $j$-th mesh node.
The set containing all our degrees-of-freedom (DoF) indices is referred to by
$\mathcal{I}_{h_{l}}$. If the multigrid level is obvious from the context, we
will try to suppress the level dependence $h_{l}$ for a more compact notation.
By $u=\sum_{i}u_{i}\phi_{i}$ and $v=\sum_{i}v_{i}\phi_{i}$ we denote linear
combinations of the nodal basis function with coefficients
$u_{i},v_{i}\in\mathbb{R}$. Defining the matrix $A_{ij}=a(\phi_{i},\phi_{j})$
and vector $\mathbf{f}_{i}=F(\phi_{i})$ results in the linear algebraic
formulation of the discrete variational problem associated with the weak
formulation: $\text{Find }\mathbf{u}\in\mathbb{R}^{|\mathcal{I}|}\text{
satisfying }A\mathbf{u}=\mathbf{f}.$

Hybrid meshes impose a domain-partitioning, which is also used for assigning
the DoF in an HPC environment to computing nodes. This approach avoids
communication between the DoF located inside the same macro-primitive, while
for DoF on different macro-primitives communication is necessary. This has to
be considered for an efficient evaluation of our operators since operations
acting locally on the same primitive type do not require inter-node
communication. To define these local operations, we have to introduce notation
to localize our vectors and matrices: For arbitrary index sets
$I\operatorname{\subseteq}\mathcal{I}$ we define restriction operators
$R_{I}:\mathbb{R}^{|\mathcal{I}|}\to\mathbb{R}^{|I|}$ consisting of zeros and
ones, which discard vector entries whose component is not present in the index
set and just retain entries in $I$. We also assume that the restriction
operator retains the global DoF ordering. Given a macro-primitive
$p\in\mathcal{P}_{H}$, we denote the set of all DoF which are located on the
primitive by $\mathcal{I}_{p}\operatorname{\subseteq}\mathcal{I}$ and its
restriction operator by $R_{p}=R_{\mathcal{I}_{p}}$. For an arbitrary macro-
tetrahedron $t\in\mathcal{T}_{H}$, which is adjacent to macro-vertices
$v_{i}\in\mathcal{V}_{H}$, $1\leq i\leq 4$, macro-edges
$e_{j}\in\mathcal{E}_{H}$, $1\leq j\leq 6$ and macro-faces
$f_{k}\in\mathcal{F}_{H}$, $1\leq k\leq 4$ we define the index-set of its
ghost-layer as $\mathcal{I}_{\partial
t}=(\cup_{i=1}^{4}\mathcal{I}_{v_{i}})\cup(\cup_{j=1}^{6}\mathcal{I}_{e_{j}})\cup(\cup_{k=1}^{4}\mathcal{I}_{f_{k}})$.
All these sets are illustrated in Figure 1.

Our surrogate ILU algorithm heavily relies on geometric properties associated
with our DoF: Each micro-vertex in a macro-tetrahedron on level $L$ can be
labeled by the logical grid coordinates
$G^{L}_{t}=\\{(x,y,z)\in\mathbb{Z}^{3}\,:\,0\leq x,y,z\textmd{ and
}x+y+z<2^{L}+1\\}$. Similarly, we define the inner grid coordinates by
$\mathring{G}^{L}_{t}=\\{(x,y,z)\in\mathbb{Z}^{3}\,:\,1\leq x,y,z\textmd{ and
}x+y+z<2^{L}\\}$. If we restrict the coordinates by setting z to a fixed
value, we obtain a face-layer $G^{N}_{f}=\\{(x,y)\in\mathbb{Z}^{2}\,:\,0\leq
x,y\textmd{ and }x+y<N\\}.$ For a vector $\mathbf{u}|_{I_{t}\cup\partial
I_{t}}$ on level $l$ restricted to a tetrahedron $t$, there is a one-to-one
correspondence between DoF-indices in $\mathcal{I}_{t}$ and inner logical grid
coordinates $\mathring{G}^{L}_{t}$ which can be constructed as follows: Assume
that $t$ is adjacent to the macro-vertices $v_{i}$ at coordinates
$\widetilde{p}_{i}\in\mathbb{R}^{3}$ for $1\leq i\leq 4$. The tetrahedron is
spanned by the edges $d_{i}=\widetilde{p}_{i+1}-\widetilde{p}_{1}$ at the base
point $\widetilde{p}_{1}$ for $1\leq i\leq 3$ (see Figure 2 left). The point
$\widetilde{p}^{(x,y,z)}=(d_{1}\cdot x+d_{2}\cdot y+d_{3}\cdot z)/(2^{L}+1)$
for $(x,y,z)\in{G}^{L}_{t}$ belongs to a shape function $\phi_{k}\in
V_{h_{l}}$ with $k\in\mathcal{I}_{t}\cup\mathcal{I}_{\partial t}$ such that
$\phi_{k}(\widetilde{p}^{(x,y,z)})=1$. This induces the mapping
$\iota_{t}:\mathcal{I}_{t}\cup\mathcal{I}_{\partial t}\to{G}^{L}_{t}$ with
$\iota_{t}(k)=(x,y,z)$. Thus, vector components $u_{i}$ of a vector
$\mathbf{u}|_{\mathcal{I}_{t}\cup\mathcal{I}_{\partial t}}$ with
$\iota\in\mathcal{I}_{t}\cup\mathcal{I}_{\partial t}$ will also be referred to
by $u^{(x,y,z)}$ or $u^{p}$ for $p=(x,y,z)\in G^{L}_{t}$, when the macro-
tetrahedron $t$ is evident from the context.

The mapping between logical grid coordinates and local DoF in
$\mathcal{I}_{t}$ also allows us to specify an ordering of the DoF indices.
This is crucial since the properties of the Gauss–Seidel-Smoother (GS-
Smoother) or the ILU-Smoother strongly depend on this order. For
$i,j\in\mathcal{I}_{t}$ with logical grid coordinates
$(x_{i},y_{i},z_{i})=\iota_{t}(i)$ and $(x_{j},y_{j},z_{j})=\iota_{t}(j),$ we
fix the ordering by

$i<j\implies(z_{i}<z_{j})\vee(z_{i}=z_{j}\wedge
y_{i}<y_{j})\vee(z_{i}=z_{j}\wedge y_{i}=y_{j}\wedge x_{i}<x_{j}).$

Consequently, the ordering strongly depends on the order of the adjacent
macro-vertices $v_{i}$ which we used to construct $\iota_{t}$. In Section 3,
we will use this by permutating $v_{i}$ with a permutation $\pi$ to obtain
good smoothing factors $\mu_{t}$ for the new DoF ordering.

Figure 2:  Left: Direction vectors in a one-to-one correspondence between DoF
and coordinates. Right: Stencil directions and grid coordinates inside a
structured tetrahedral grid.

We will now introduce the applied stencil notation for our surrogate-ILU-
Algorithm in Section 4. For this, we first define the stencil directions
between logical coordinates as displacement vectors, i.e. $\\{x-y\,|\,x,y\in
G^{L}_{t}\\}$. The most common directions are named after the four cardinal
directions, as well as the top and bottom directions, such that the x-axis
runs from west to east, the y-axis from south to north, and the z-axis from
top to bottom. For instance, the west direction $w$ corresponds to the
displacement $(-1,0,0)$. All stencil directions are collected in the set

$\mathcal{D}=\\{w,s,se,bnw,bn,bc,be,c,e,n,nw,tse,ts,tc,tw\\}.$

Relying on the ordering defined above, the set of all lower stencil directions
needed in our ILU is given by

$\mathcal{D}_{l}=\\{w,s,se,bnw,bn,bc,be\\}.$

Consider two indices $i\in\mathcal{I}_{t}$ and
$j\in\mathcal{I}_{t}\cup\mathcal{I}_{\partial t}$ with coordinates
$p_{i}=\iota_{t}(i)$ and $p_{j}=\iota_{t}(j)$. Due to the local support of the
low order conforming finite-element shape functions, we know that if
$A_{ij}\neq 0$ there exists a $\widetilde{d}\in\mathcal{D}$ such that
$p_{j}=p_{i}+\widetilde{d}$. We can define the stencil
$(A^{p_{i}}_{d})_{d\in\mathcal{D}}$ by $A^{p_{i}}_{\widetilde{d}}=A_{ij}$. The
matrix-vector multiplication $v|_{\mathcal{I}_{t}}=(Au)|_{\mathcal{I}_{t}}$ on
the macro-tetrahedron $t$ can therefore be written in terms of stencils as

$v^{p}=\sum_{d\in\mathcal{D}}A^{p}_{d}u^{p+d}\quad\textmd{ for all }\quad
p\in\mathring{G}^{L}_{t},$

where we identified the DoFs with logical coordinates. Stencils and the
associated grid coordinates are depicted in Figure 2 (right).

We mainly rely on a geometric multigrid algorithm which combines a so called
smoother with a coarse grid correction step to an optimal solver (see
e.g.[23]). We now introduce smoothers which only act on the DoF of a single
primitive. This is motivated by our hybrid mesh on which only operations
between DoF located on the same primitive are cheap while everything else
requires expensive inter-node communication. In our setting, for a given
primitive $p\in\mathcal{P}_{H}$, a smoother acting on the DoF
$\mathcal{I}_{p}$ located on the primitive can be described by applying a
preconditioner matrix
$C_{p}\in\mathbb{R}^{|\mathcal{I}_{p}|\times|\mathcal{I}_{p}|}$ inside a
Richardson iteration with the appropriate restriction operators

$\mathbf{u}\leftarrow\mathbf{u}+R^{T}_{p}C^{-1}_{p}R_{p}(\mathbf{f}-A\mathbf{u}),$

where we denote the current right-hand-side vector by $\mathbf{f}$ and the
current estimate by $\mathbf{u}$. Note that usually $C_{p}$ depends on our
system matrix $A$ and is not necessarily symmetric. If it is, we stress this
by using $C_{p,sym}$.

As a reference example, we consider the GS-Smoother: On the primitive $p$ we
define the lower triangular part
$L_{A,p}\in\mathbb{R}^{|\mathcal{I}_{p}|\times|\mathcal{I}_{p}|}$ and diagonal
part $D_{A,p}\in\mathbb{R}^{|\mathcal{I}_{p}|\times|\mathcal{I}_{p}|}$ of our
system matrix $A$ restricted to the primitive as

$(L_{A,p})_{ij}=\begin{cases}(R_{p}AR^{T}_{p})_{ij}&\textmd{ if }j\leq i\\\
0&\textmd{ else }\end{cases}\quad\textmd{ and
}\quad(D_{A,p})_{ij}=\begin{cases}(R_{p}AR^{T}_{p})_{ii}&\textmd{ if }i=j\\\
0&\textmd{ else }\end{cases}.$

The GS-Smoother is then given by $C_{p}=L_{A,p}$ and its symmetrized version
by $C_{p,sym}=(L_{A,p})(D_{A,p})^{-1}(L_{A,p})^{T}$.

###  2.2 Steklov–Poincaré operator

The idea of many non-overlapping domain decomposition methods is to solve a
system of equations just on the boundary primitives of a decomposition and
apply static condensation to the interior DoF. Often, the decomposition is
performed on unstructured meshes using graph partitioning libraries like METIS
[25] or SCOTCH [36]. On hybrid meshes, the decomposition can simply be derived
from the macro-mesh such that each refined macro-tetrahedron represents a
subdomain. The macro-vertices, edges and faces and their respective DoF then
form the interfaces between the subdomains. The index set containing all DoF
on the macro-interface is given by
$\mathcal{I}_{\Gamma}=\cup_{t\in\mathcal{T}_{H}}\mathcal{I}_{\partial t}$.
With the previously introduced restriction operators, we define the submatrix
$A_{\Gamma\Gamma}=R_{\Gamma}AR^{T}_{\Gamma}$ which just considers the coupling
between the boundary DoFs, the matrix $A_{tt}=R_{t}AR^{T}_{t}$ which couples
the interior DoF of a macro-tetrahedron $t\in\mathcal{T}_{h}$ and the coupling
matrices $A_{\Gamma t}=R_{\Gamma}AR^{T}_{t}$,
$A_{t\Gamma}=R_{t}AR^{T}_{\Gamma}$ between the interior DoF of a tetrahedron
and the boundary DoF. Inverting $A$ is equivalent to solving
$S\mathbf{u}_{\Gamma}=\bm{\chi}_{\Gamma}$ for $\mathbf{u}_{\Gamma}$ where the
Steklov–Poincaré operator $S$ is given by

$S=A_{\Gamma\Gamma}-\sum_{t\in\mathcal{T}_{H}}A_{\Gamma{}t}A_{tt}^{-1}A_{t\Gamma}$

and the right hand side is
$\bm{\chi}_{\Gamma}=R_{\Gamma}\mathbf{b}-\sum_{t\in\mathcal{T}_{H}}A_{\Gamma
t}A^{-1}_{tt}R_{t}\mathbf{b}$. The interior DoFs can be reconstructed by
solving $A_{tt}\mathbf{u}_{t}=R_{t}\mathbf{b}-A_{t\Gamma}\mathbf{u}_{\Gamma}$
for $\mathbf{u}_{t}$ on each tetrahedron $t\in\mathcal{T}_{H}$. Usually, $S$
is not constructed explicitly, but is inverted by using a preconditioned
conjugate gradient (PCG) method (see [15, 43] for a non-exhaustive overview).
The PCG algorithms relies on the evaluation of matrix-vector products with
$S$. Evaluating $S$ requires the evaluation of $A_{tt}^{-1}$ for which we use
a multigrid method with an ILU-Smoother in the interior of our macro-
tetrahedron. We will show, that our ILU algorithm gives us robustness with
respect to tetrahedra which are distorted along one axial direction.

##  3 Matrix-based ILU-Smoother

In this section, we introduce an ILU-Smoother that is amenable to an efficient
matrix-free implementation and investigate its performance. Both the
efficiency and the possibility for a matrix-free implementation heavily depend
on the used ILU formulation.

We first start with the general definition: Let the sparsity pattern of the
matrix $A$ be given by
$\mathcal{S}_{A}=\\{(i,j)\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{I}|}|A_{ij}\neq
0\\},$ then the ILU factorization consisting of the lower-triangular matrix
$L\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{I}|}$ and diagonal matrix
$D\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{I}|}$ is defined by
$(LDL^{T})_{ij}=A_{ij}$ for $i,j\in\mathcal{S}_{A}$.

The ILU will be restricted to the interior of macro-tetrahedra. Hence, we
define the lower triangular matrices
$L_{t}\in\mathbb{R}^{|\mathcal{I}_{t}|\times|\mathcal{I}_{t}|}$ and diagonal
matrices $D_{t}\in\mathbb{R}^{|\mathcal{I}_{t}|\times|\mathcal{I}_{t}|}$ by
$(L_{t}D_{t}L_{t}^{T})_{ij}=(R_{t}AR^{T}_{t})_{ij}$ for
$i,j\in\mathcal{S}_{R_{t}AR^{T}_{t}}$. We thus define our preconditioner on
the macro-tetrahedra by $C_{t,sym}=L_{t}D_{t}L^{T}_{t}$.

######  Remark 3.1.

Another popular choice would be the modified-ILU which is derived from the
factorization $(\mathring{L}+D)D^{-1}(\mathring{L}+D)^{T}_{ij}=A_{ij}$ for
$(i,j)\in\mathcal{S}_{A}$, where $\mathring{L}$ is now a strictly lower-
diagonal matrix. The remainder matrix containing the additional fill-in is
$R=(\mathring{L}+D)D^{-1}(\mathring{L}+D)-A$. Summing over the remainder
matrix terms and adding them to the diagonal with a weight $\omega$, i.e.
$(D_{\omega})_{ii}=D_{ii}+\omega\sum_{i\neq j}\left|r_{ij}\right|$, yields the
modified ILU

$(\mathring{L}+D_{\omega})D_{\omega}^{-1}(\mathring{L}+D_{\omega})^{T},$

which for $\omega=0$ becomes the usual ILU factorization, but is known to
behave more robustly for different triangle types [35]. For $\omega>0$, the
$D_{\omega}$ matrix contains entries from the fill-in compared to the sparse
matrix $A$, which we aim to approximate with the ILU. In 3D, our fill-in
consists of 12 additional nonzero entries per row which cannot be calculated
in a memory-efficient way during our factorization. They can possibly be
reconstructed approximately in a postprocessing step, but we nevertheless
restrict ourselves to the case of $\omega=0$.

For our concrete implementation, we aim for a factorization of the form
$LDL^{T}$, where $L$ is a lower triangular matrix with a unit diagonal to
minimize the number of multiplications and divisions during the forward and
backward substitutions.

###  3.1 Strategy on macro-elements

The performance of the ILU smoother strongly depends on the ordering of the
DoF. To mitigate this effect and increase robustness, an alternating ILU-
Smoother can be used [35] which successively applies several ILU
factorizations with different orderings. This cannot be efficiently done in a
matrix-free algorithm since the different orderings would result in cache
misses for at least one of the orderings. Thus, we have to use an efficient
ordering from the beginning which means that we may have to permute the
tetrahedral vertices in a preprocessing step. It is shown in [37] with a Local
Fourier Analysis (LFA) that for optimal performance, the triangles in 2D have
to be orientated such that the first vertex is at the largest angle, the
second at the smallest, and the third at the remaining angle.

This heuristic does not directly extend to 3D. For a more reliable strategy,
we iterate over all macro-tetrahedra and apply an LFA with the techniques from
[20, 21] to the asymptotic ILU stencil for each of its possible orientations.
Finally, we apply the permutation resulting in the smallest smoothing factor
to the macro-tetrahedron.

For a formula of the asymptotic ILU stencil, we derive the ILU in-place
factorization which only relies on local information from DoF neighbors. On
structured grids, the factorization can be derived by evaluating $LDL^{T}$ and
comparing it to the matrix $A$ on the respective sparsity pattern. This yields
the stencil equations

$A^{p}_{d}=\sum_{\widetilde{d}\in\mathcal{D}_{l}\cup\\{c\\}}L^{p}_{\widetilde{d}}D^{p+\widetilde{d}}_{c}L^{p+d}_{\widetilde{d}-d},\quad\textmd{with}\quad
d\in\mathcal{D}_{l}\cup\\{c\\}.$

If the stencils at $p+d$ for $d\in\mathcal{D}_{l}$ are already factorized,
this system of equations can be solved for $L^{p}_{\widetilde{d}}$ and
$D^{p}_{c}$ for $\widetilde{d}\in\mathcal{D}_{l}\cup\\{c\\}$. Thus, the
factorization of the stencils $L^{p}_{\widetilde{d}}$ and $D^{p}_{c}$ at
$p=(x,y,z)$ at level $L$ relies on information at the points

(2) $\displaystyle I_{p}$ $\displaystyle=I^{\beta}_{p}\cup
I^{\gamma}_{p}\quad\textmd{where}$ (3) $\displaystyle I^{\beta}_{p}$
$\displaystyle=\left\\{(\widetilde{x},y,z)|0\leq\widetilde{x}<x\right\\}\cup\left\\{(\widetilde{x},\widetilde{y},z)|0\leq\widetilde{x}<2^{L}\\!+\\!1\\!-\\!\widetilde{y}\\!-\\!z,\widetilde{y}<y\right\\}\quad\textmd{
and }$ (4) $\displaystyle I^{\gamma}_{p}$
$\displaystyle=\left\\{(\widetilde{x},\widetilde{y},z-1)|0\leq\widetilde{x}<2^{L}\\!+\\!2\\!-\\!\widetilde{y}\\!-\\!z,0\leq\widetilde{y}<2^{L}\\!+\\!2\\!-\\!z\right\\}$

are the logical DoF-coordinates on the current layer and the layer below. This
standard procedure was already applied in the very first ILU paper[34] for an
incomplete Cholesky decomposition in 2D on quadrilateral grids.

Specifically, we obtain the following equations

(5) $\displaystyle A^{p}_{bc}$ $\displaystyle=L^{p}_{bc}\bm{D^{p+bc}_{c}},$
(6) $\displaystyle A^{p}_{s}$
$\displaystyle=L^{p}_{bc}\bm{D^{p+bc}_{c}}\bm{L^{p+s}_{bn}}+L^{p}_{s}\bm{D^{p+s}_{c}},$
(7) $\displaystyle A^{p}_{bnw}$
$\displaystyle=L^{p}_{bc}\bm{D^{p+bc}_{c}}\bm{L^{p+bnw}_{se}}+L^{p}_{bnw}\bm{D^{p+bnw}_{c}},$
(8) $\displaystyle A^{p}_{be}$
$\displaystyle=L^{p}_{bc}\bm{D^{p+bc}_{c}}\bm{L^{p+be}_{w}}+L^{p}_{be}\bm{D^{p+be}_{c}},$
(9) $\displaystyle A^{p}_{w}$
$\displaystyle=L^{p}_{bc}\bm{D^{p+bc}_{c}}\bm{L^{p+w}_{be}}+L^{p}_{bnw}\bm{D^{p+bnw}_{c}}\bm{L^{p+w}_{bn}}+L^{p}_{s}\bm{D^{p+s}_{c}}\bm{L^{p+w}_{se}}+L^{p}_{w}\bm{D^{p+w}_{c}},$
(10) $\displaystyle A^{p}_{bn}$
$\displaystyle=L^{p}_{bc}\bm{D^{p+bc}_{c}}\bm{L^{p+bn}_{s}}+L^{p}_{be}\bm{D^{p+be}_{c}}\bm{L^{p+bn}_{se}}+L^{p}_{bnw}\bm{D^{p+bnw}_{c}}\bm{L^{p+bn}_{w}}+L^{p}_{bn}\bm{D^{p+bn}_{c}},$
(11) $\displaystyle A^{p}_{se}$
$\displaystyle=L^{p}_{bc}\bm{D^{p+bc}_{c}}\bm{L^{p+se}_{bnw}}+L^{p}_{be}\bm{D^{p+be}_{c}}\bm{L^{p+se}_{bn}}+L^{p}_{se}\bm{D^{p+se}_{c}}+L^{p}_{s}\bm{D^{p+s}_{c}}\bm{L^{p+se}_{w}}$
and (12) $\displaystyle A^{p}_{c}$
$\displaystyle=D^{p}_{c}+\left(L^{p}_{bc}\right)^{2}\bm{D^{p+bc}_{c}}+\left(L^{p}_{be}\right)^{2}\bm{D^{p+be}_{c}}+\left(L^{p}_{bnw}\right)^{2}\bm{D^{p+bnw}_{c}}+\left(L^{p}_{bn}\right)^{2}\bm{D^{p+bn}_{c}}$
(13)
$\displaystyle\qquad+\left(L^{p}_{se}\right)^{2}\bm{D^{p+se}_{c}}+\left(L^{p}_{s}\right)^{2}\bm{D^{p+s}_{c}}+\left(L^{p}_{w}\right)^{2}\bm{D^{p+w}_{c}},$

at a grid point $p=(x,y,z)$, where the already factorized symbols at points in
$I_{p}$ are formatted in bold for convenience. We further omitted the
equations which require information that is not available during the
factorization.

To estimate the asymptotic stencils, we assume $K=\textmd{Id}_{3}$ inside the
macro-tetrahedra. The assumption is valid if $K$ is close to a scaled identity
and it is not spatially varying too much. Inside the macro-tetrahedra, this
leads to coordinate independent stencils $A^{p}_{d}=A_{d}$ for all
$d\in\mathcal{D}$ at all grid points $p$. Similarly to the 2D case from [37],
we can now calculate the asymptotic ILU stencils $L^{\infty}_{d}$ and
$D^{\infty}_{c}$ for $d\in\mathcal{D}$ by iterating Eqs. 13 with a Gauss-
Seidel scheme. The initial stencil values are given by
$L^{(0)}_{c}=D^{(0)}_{c}=1,L^{(0)}_{d}=D^{(0)}_{d}=0$ for
$d\in\mathcal{D}\setminus\\{c\\}$.

The Fourier symbols are then given by

(14) $\displaystyle D({\bm{\theta}})=D^{\infty}_{c},\quad
L({\bm{\theta}})=1+\sum_{d\in\mathcal{D}_{l}}L^{\infty}_{d}e^{id\cdot{\bm{\theta}}}\quad\textmd{
and }\quad
A({\bm{\theta}})=\sum_{d\in\mathcal{D}}A_{d}e^{id\cdot{\bm{\theta}}}$

for ${\bm{\theta}}\in\left(-\pi,+\pi\right)^{3}$. The smoothing factor
$\mu_{t}$ is then calculated for the highly-oscillating frequencies
$\Theta_{osc}=\left(-\pi,+\pi\right)^{3}\setminus\left(-\pi/2,+\pi/2\right)^{3}$
by

$\mu_{t}=\sup_{{\bm{\theta}}\in\Theta_{osc}}\left|\frac{L({\bm{\theta}})D({\bm{\theta}})\conj{L}({\bm{\theta}})-A({\bm{\theta}})}{L({\bm{\theta}})D({\bm{\theta}})\conj{L}({\bm{\theta}})}\right|.$

For our purpose, $\mu_{t}$ is determined by sampling $\Theta_{osc}$ on a
uniform grid with $16$ samples in each coordinate direction.

In order to determine the optimal orientation of a macro-tetrahedron $t$, we
iterate over all possible vertex permutations $\pi_{k}$, $1\leq k\leq 24$,
calculate $\mu_{t_{\pi_{k}}}$ for all the permutated tetrahedra $t_{\pi_{k}}$
and finally apply the permutation with the smallest smoothing factor
$\mu_{t_{\pi_{k}}}$. Note that we need to calculate $16^{3}\cdot 24=98304$
symbols per macro-tetrahedron to determine the final orientation, but we only
have to do this once during a preprocessing step before the actual
calculation.

######  Remark 3.2.

We can a-priori determine the scaling behaviour of the asymptotic stencils in
the mesh width $h$: For this, assume that the asymptotic stencils
$D^{\infty}_{c}$, $L^{\infty}_{d}$ for $d\in\mathcal{D}_{l}$ exist and that
the stencil $A^{p}_{d}$ is independent of $p$. Due to the transformation rule,
we know $A_{d}\sim\Theta(h)$ and since
$A_{c}=D^{\infty}_{c}(1+\sum_{d\in\mathcal{D}_{l}}(L^{\infty}_{d})^{2})$ has
just positive summands, we immediately see $D^{\infty}_{c}\sim\Theta(h)$,
$D^{\infty}_{c}(L^{\infty}_{d})^{2}\sim\Theta(h)$ and thus
$L^{\infty}_{d}\sim\Theta(1)$.

Numerical results for the performance of our smoother for different
permutations validating this strategy are given in Appendix A.

##  4 Matrix-free ILU-Smoother on macro-tetrahedra

After introducing our matrix-based hybrid ILU-Smoother, we now turn towards
our matrix-free algorithm. Our goal is to develop a surrogate smoother with
the same smoothing performance as the ILU, but with significantly lower memory
requirements supporting large scale computations. For this purpose, it is
sufficient to consider structured meshes consisting of a single macro-
tetrahedron since its behavior on hybrid grids is the same as for the matrix-
based variant. We also reflect this in our notation by omitting references to
specific tetrahedra $t\in\mathcal{T}_{H}$ whenever possible and by keeping it
as concise as possible. Furthermore, we always assume that this macro-
tetrahedron is oriented with respect to our reordering strategy.

Figure 3:  Stencil functions for the $L^{(\cdot)}_{d}$ and $D^{(\cdot)}_{c}$
stencils of the distorted tetrahedron of height $h=0.1$ on grid level 5. The
tetrahedron is scaled in the z-direction for better visibility. The original
tetrahedron is depicted with red dashed lines in the plots.

To motivate our approach, we depict in Figure 3 each stencil direction of the
factorized ILU of the distorted tetrahedron of height $h=0.1$ as continuous
functions. The z-axis is exaggerated by a factor of 10 to make the stencil
plots easier to read. The outlines of the distorted untransformed tetrahedron
are indicated by red dashed lines in the plots while the outlines of the
exaggerated tetrahedron are indicated by black lines. One vertical and three
horizontal slices make the inside visible. For all stencil directions, we see
that the stencil function does not vary much in the x-y-plane, except for some
tiny layer close to the boundary. Along the z-axis, all stencil functions
display a color gradient from top to bottom.

This allows us to make several important observations: Firstly, the stencil
weights of a single direction are smooth functions. Secondly, the stencil
weights are anisotropic with respect to the z-axis. Both observations will
enter into our surrogate strategy.

Figure 4:  Stencil plots for the distorted tetrahedron on different grid
levels along the line parallel to the z-axis starting at $(0.1,0.1,0)$ and
ending at $(0.1,0.1,0.1)$.

In Figure 4, we have a more thorough look at the behavior in the z-direction.
We have depicted the stencil weights of the matrix-based ILU for the distorted
tetrahedron along a line starting at $(0.1,0.1,0)$ and ending at
$(0.1,0.1,0.1)$ for refinement levels from $5$ to $7$. Each stencil direction
approaches its asymptotic value towards the right. From Remark 3.2, we expect
that the diagonal stencil values from the $D$ matrix scale as the diagonal
entries of the Laplacian with respect to the mesh width $h$. Thus, its
asymptotic value from a coarse level to its refinement has to decrease by a
factor of $\nicefrac{{1}}{{2}}$ which is visible in the y-scaling of the
center plots. For the stencils of the $L$ matrix, Remark 3.2 suggests that
they do not scale in $h$ at all, which results in the same lower and upper
bounds in the plots. We see that all the curves for the different levels
retain the same shape. Furthermore, the curve for level 5 scaled by
$\nicefrac{{1}}{{2}}$ or $\nicefrac{{1}}{{4}}$ in the x-direction yields
curves similar to the ones on levels 6 and 7. This suggests, that the non-
asymptotic part of the stencil factorization decays in a level-independent way
and only depends on the graph distance on the mesh to the DoF in the x-y-
plane.

###  4.1 The surrogates

We want to replace our in-memory ILU algorithm with a matrix-free variant to
save main-memory storage capacities and enable large-scale computations. Note
that we cannot use an on-the-fly approach to execute the ILU Algorithm, since
the factorization itself does not rely only on local information. Furthermore,
only the forward substitution moves through the DoF in the same direction as
the typical in-place ILU factorization algorithm while the backward
substitution moves in the opposite direction. Rewriting the factorization such
that it moves in the same direction as the backward substitution is possible,
but it is not numerically stable due to the fact that it requires to
reconstruct ILU stencils from nearly asymptotic ILU stencils. Therefore, a
memory-efficient factorization will only be the first step of our algorithm
and has to be combined with some approximation approach that reconstructs the
ILU.

For the factorization, we will use an in-place approach and store just enough
information to complete the factorization. Matrix entries which are no longer
needed for the factorization are immediately discarded. These types of in-
place factorizations are typically used for ILU factorizations[34] and just
rely on local information from factorized neighboring DoF. The factorization
can be derived by solving Eqs 13 for $L^{p}_{d},d\in\mathcal{D}_{l}$ and
$D^{p}_{c}$ at each point $p=(x,y,z)$. Our asymptotic storage requirements for
the factorization can be estimated from the set $I_{p}$. Clearly, our memory
requirements scale as $|I_{p}|=\mathcal{O}(h^{-2})$ for storing the
factorization on the current face-layer ($I_{\beta}$) and the layer below
($I_{\gamma}$). This has to be compared to a memory consumption scaling as
$\mathcal{O}(h^{-3})$ which we would need to store the stencils on the whole
grid in a matrix-based implementation. The $\mathcal{O}(h^{-2})$ for the
factorization are thereby negligible for practical purposes.

We want to replace factorized stencils with surrogate polynomial
approximations in the anisotropic polynomial space

$\mathbb{P}_{\textmd{dg}_{x},\textmd{dg}_{y},\textmd{dg}_{z}}=\left\\{x^{i}y^{j}z^{k}\,|\,0\leq
i\leq\textmd{dg}_{x},0\leq j\leq\textmd{dg}_{y},0\leq
k\leq\textmd{dg}_{z}\right\\}\,\textmd{for}\,\textmd{dg}_{x},\textmd{dg}_{y},\textmd{dg}_{z}\in\mathbb{N}$

such that
$L_{d}^{(\cdot)},(D^{-1}_{c})^{(\cdot)}\in\mathbb{P}_{\textmd{dg}_{x},\textmd{dg}_{y},\textmd{dg}_{z}}$
for $d\in\mathcal{D}_{l}\cup\\{c\\}$. The different degrees for the different
directions will allow us to vary the accuracy of the stencils in a direction-
dependent way. The factorized stencils on our grid will be added direction-
wise to a least-squares (LSQ) problem to obtain one surrogate polynomial per
stencil direction. Of course, adding all the stencils to our least-squares
problem would require a large memory overhead to store the matrix and a large
computational overhead for solving the least-squares problem. Therefore, only
stencils $L^{p}_{d}$, $D^{p}_{c}$ at points $p\in S_{d}$ on some coarser grid
$S_{d}$ of sample points will be added whose structure depends on the
direction $d\in\mathcal{D}_{l}\cup\\{c\\}$. For a coarse sampling grid at
level $L_{H}$, these sample points are defined as

$S_{d}=\left\\{p\,|\,p\in
G^{L}_{t},(p-(\mathbf{1}-d))\,\textmd{mod}\,\operatorname*{\vphantom{p}min}(2^{L}/2^{L_{H}},1)=0\textmd{
and }p+d\not\in G^{L}_{t}\setminus\mathring{G}^{L}_{t}\right\\}.$

They are illustrated for three directions in Figure 5.

Figure 5:  Illustration of the sample sets in the x-z plane for $be$, $w$, and
$c$. All points have the distance $2^{L}/2^{L_{H}}$ along each grid axis. Note
the offset from the bottom for $be$ and the offset from the left for $w$. No
shift has to be applied for the $c$ direction.

The shift by direction $d$ is necessary to make sure to have the utmost
boundary points in our approximation as sample points, which impacts the
quality of the surrogate and thus the smoother.

1:$N=\textmd{ number of micro-vertices on a macro edge}$.

2:Initialize $\beta^{(x,y)}_{d}=\gamma^{(x,y)}_{d}=0$ for
$d\in\mathcal{D}_{l}$ and $\beta^{(x,y)}_{c}=\gamma^{(x,y)}_{c}=1$ and
$(x,y)\in G_{f}^{N}$.

3:for $z=1,\dots,N-2$ do

4: for $y=1,\dots,N-2-z$ do

5: for $x=1,\dots,N-2-z-y$ do

6: Assemble $A_{d}=A^{(x,y,z)}_{d}$ for $d\in\mathcal{D}$.

7: Assign

8: $\begin{aligned} \beta^{(x,y)}_{bc}&=A_{bc}/\bm{\upgamma^{(x,y)}_{c}}\\\
\beta^{(x,y)}_{s}&=(A_{s}-{\beta^{(x,y)}_{bc}}\bm{\upgamma^{(x,y)}_{c}}\bm{\upbeta^{(x,y-1)}_{bn}})/\bm{\upbeta^{(x,y-1)}_{c}}\\\
\beta^{(x,y)}_{bnw}&=(A_{bnw}-{\beta^{(x,y)}_{bc}}\bm{\upgamma^{(x,y)}_{c}}\bm{\upgamma^{(x-1,y+1)}_{se}})/\bm{\upgamma^{(x-1,y+1)}_{c}}\\\
\beta^{(x,y)}_{be}&=(A_{be}-{\beta^{(x,y)}_{bc}}\bm{\upgamma^{(x,y)}_{c}}\bm{\upgamma^{(x+1,y)}_{w}})/\bm{\upgamma^{(x+1,y)}_{c}}\\\
\beta^{(x,y)}_{w}&=(A_{w}-{\beta^{(x,y)}_{bc}}\bm{\upgamma^{(x,y)}_{c}}\bm{\upbeta^{(x-1,y)}_{be}}-{\beta^{(x,y)}_{bnw}}\bm{\upgamma^{(x-1,y+1)}_{c}}\bm{\upbeta^{(x-1,y)}_{bn}}\\\
&\qquad-{\beta^{(x,y)}_{s}}\bm{\upbeta^{(x,y-1)}_{c}}\bm{\upbeta^{(x-1,y)}_{se}})/\bm{\upbeta^{(x-1,y)}_{c}}\\\
\beta^{(x,y)}_{bn}&=(A_{bn}-{\beta^{(x,y)}_{bc}}\bm{\upgamma^{(x,y)}_{c}}\bm{\upgamma^{(x,y+1)}_{s}}-{\beta^{(x,y)}_{be}}\bm{\upgamma^{(x+1,y)}_{c}}\bm{\upgamma^{(x,y+1)}_{se}}\\\
&\qquad-{\beta^{(x,y)}_{bnw}}\bm{\upgamma^{(x-1,y+1)}_{c}}\bm{\upgamma^{(x,y+1)}_{w}})/\bm{\upgamma^{(x,y+1)}_{c}}\\\
\beta^{(x,y)}_{se}&=(A_{se}-{\beta^{(x,y)}_{bc}}\bm{\upgamma^{(x,y)}_{c}}\bm{\upbeta^{(x+1,y-1)}_{bnw}}-{\beta^{(x,y)}_{be}}\bm{\upgamma^{(x+1,y)}_{c}}\bm{\upbeta^{(x+1,y-1)}_{bn}}\\\
&\qquad-{\beta^{(x,y)}_{s}}\bm{\upbeta^{(x,y-1)}_{c}}\bm{\upbeta^{(x+1,y-1)}_{w}})/\bm{\upbeta^{(x+1,y-1)}_{c}}\\\
\beta^{(x,y)}_{c}&=A_{c}-\left({\beta^{(x,y)}_{bc}}\right)^{2}\bm{\upgamma^{(x,y)}_{c}}-\left({\beta^{(x,y)}_{be}}\right)^{2}\bm{\upgamma^{(x+1,y)}_{c}}\\\
&\qquad-\left({\beta^{(x,y)}_{bnw}}\right)^{2}\bm{\upgamma^{(x-1,y+1)}_{c}}-\left({\beta^{(x,y)}_{bn}}\right)^{2}\bm{\upgamma^{(x,y+1)}_{c}}\\\
&\qquad-\left({\beta^{(x,y)}_{se}}\right)^{2}\bm{\upbeta^{(x+1,y-1)}_{c}}-\left({\beta^{(x,y)}_{s}}\right)^{2}\bm{\upbeta^{(x,y-1)}_{c}}\\\
&\qquad-\left({\beta^{(x,y)}_{w}}\right)^{2}\bm{\upbeta^{(x-1,y)}_{c}}\end{aligned}$

9: Set $p=h\cdot(x,y,z)\in\mathbb{R}^{3}$.

10: for $d\in\mathcal{D}_{l}$ do

11: if $(x,y,z)\in S_{d}$ then Add $\beta^{(x,y)}_{d}$ at $p$ to a LSQ
problem.

12: if $(x,y,z)$ on the cell boundary then Store $\beta^{(x,y)}_{d}$ in
memory.$\triangleright$ (V1)

13: if $(x,y,z)\in S_{c}$ then Add $1/\beta^{(x,y)}_{c}$ at $p$ to a LSQ
problem.

14: if $(x,y,z)$ on the cell boundary then Store $1/\beta^{(x,y)}_{c}$ in
memory.$\triangleright$ (V1)

15: Copy $\beta$ into $\gamma$.

16: Set $\beta^{(x,y)}_{c}=1$ and $\beta^{(x,y)}_{d}=0$ for
$d\in\mathcal{D}_{l}$ with $(x,y)\in G_{f}^{N-2-z}$.

17:Solve the least-squares problems.

Algorithm 1 Memory efficient in-place $LDL^{T}$ factorization in a single
tetrahedron.

This discontinuity also affects our substitution at the boundary. For these
stencils, we will test two variants: In the first variant denoted by (V1), all
the boundary stencils $L^{p}_{d},D^{p}_{d}$ at boundary points $p$ will be
stored and later used in the backward and forward substitutions. Again, this
just requires a memory complexity of $\mathcal{O}(h^{-2})$ since we have four
boundary faces each with a $\mathcal{O}(h^{-2})$ storage requirement. In a
second variant (V2), we will refrain from storing any boundary stencils and
just correct the interpolated stencil entries at the boundary during the
substitutions by setting the respective directions to 0.

The full factorization algorithm for 3D is given in Alg. 1. All the
information in $I_{p}$ from previously factorized stencils is colored in red.
The information of the current stencil, which was calculated before and is
used for other stencil entries, is colored in blue. We denote temporarily
saved stencils of the current face-layer belonging to $I^{\beta}$ by $\beta$
and stencils on the previous face-layer which correspond to $I^{\gamma}$ by
$\gamma$. After each factorization step, we check if the stencil direction
should be added to a least-squares problem. Note that for the central
direction, we add the multiplicative inverse to the least-squares problem to
avoid a floating-point division while applying the smoother.

Note that we rescale the integer coordinates $(x,y,z)$ by the mesh width $h$
before adding them to the least-squares problem in order to have a level
independent approximation quality. The steps in which we store the boundary
stencils of the factorization are marked by (V1) and they are only used in the
first variant of the algorithm. Before increasing the z index, we have to copy
$\beta$ into $\gamma$.

1:for $z=1,\dots N-2$ do

2: for $y=1,\dots N-2-z$ do

3: for $x=1,\dots N-2-z-y$ do

4: Set $p=(x,y,z)$.

5: Evaluate $A^{p}_{d}$ for $d\in\mathcal{D}$ by surrogates with Newton’s
Divided Differences

6: Formula (NDDF) or assemble it.

7: Apply $w^{p}=b^{p}-\sum_{d\in\mathcal{D}}A^{p}_{d}x^{p+d}$ for
$d\in\mathcal{D}_{l}$.

8: if $p$ is on the cell boundary then

9: Load $L^{p}_{d}$ for $d\in\mathcal{D}_{l}$ from main memory.
$\triangleright$ (V1)

10: If $p+d$ is on the boundary set $L^{p}_{d}=0$. $\triangleright$ (V2)

11: else

12: Evaluate $L^{p}_{d}$ for $d\in\mathcal{D}_{l}$ by surrogates with NDDF.

13: $w^{p}=w^{p}-\sum_{d\in\mathcal{D}_{l}}L^{p}_{d}w^{p+d}$

14:for $z=N-2,\dots 1$ do

15: for $y=N-2-z,\dots 1$ do

16: for $x=N-2-z-y,\dots 1,$ do

17: Set $p=(x,y,z)$

18: Evaluate $(D^{p}_{c})^{-1}$ by surrogates with NDDF.

19: $w^{p}=w^{p}\cdot(D^{p}_{c})^{-1}$

20: for $d\in\mathcal{D}_{l}$ do

21: if $p-d$ is on the cell boundary then

22: Load $L^{p-d}_{d}$ for $d\in\mathcal{D}_{l}$ from main memory.
$\triangleright$ (V1)

23: If $p-d$ is on the boundary set $L^{p-d}_{d}=0$ for $d\in\mathcal{D}_{l}$.
$\triangleright$ (V2)

24: else

25: Evaluate $L^{p-d}_{d}$ for $d\in\mathcal{D}_{l}$ by surrogates with NDDF.

26: $w^{p}=w^{p}-\sum_{d\in\mathcal{D}_{l}}L^{p-d}_{d}w^{p-d}$

27: $x^{p}=x^{p}+w^{p}$

Algorithm 2 Surrogate evaluation for $x+(LDL^{T})^{-1}(b-Ax)$ with boundary
layer.

In Alg. 2, we show one step of the smoothing algorithm. Again, the steps
specific to the first variant of the algorithm are marked by (V1) while steps
exclusively for the second variant are marked by (V2).

Altogether, by merging the operator application with the forward substitution
and the backward substitution with the diagonal scaling, we have to iterate
twice over the entire mesh. Therefore, the runtime due to the storage accesses
should be comparable to the symmetric GS algorithm which also iterates twice
over the tetrahedral grid. During the first iteration, we calculate the
residual by either assembling the $A$ matrix or evaluating its surrogate
matrix. We then plug the residual into the forward substitution for $L$. In
the second iteration, we combine the multiplication of the inverse diagonal
matrix $D^{-1}$ with the backward substitution of $L^{T}$ and the correction
step. The evaluation of all the surrogate polynomials is implemented with a
Newton’s Divided Differences Formula (NDDF). The algorithm is known to be
unstable, however, this does not pose a problem for the low polynomial degrees
considered here. This approach significantly speeds up the polynomial
evaluation compared to other schemes: If the stencils for $L$ at $p=(x,y,z)$
were already evaluated, a new evaluation at the direct neighbor $(x+1,y,z)$
only takes
$(\textmd{dg}_{x}+1)\cdot|\mathcal{D}_{l}|=7\cdot(\textmd{dg}_{x}+1)$ floating
point additions and no multiplication at all. Similarly, if during the second
iteration the stencils for $L^{T}$ and $D^{-1}$ at $p=(x,y,z)$ were already
evaluated, a new evaluation at the direct neighbor $(x+1,y,z)$ only takes
$(\textmd{dg}_{x}+1)\cdot(|\mathcal{D}_{l}|+1)=8\cdot(\textmd{dg}_{x}+1)$
floating point additions. Since the number of stencils with a predecessor to
the west is $\mathcal{O}(h^{-3})$ compared to $\mathcal{O}(h^{-2})$ stencils
without a predecessor, the surrogate evaluation asymptotically becomes
computationally cheap for small $h$ and large multigrid levels. Furthermore,
the stencils $L^{p}$, $D^{p}$ and $(L^{T})^{p}$ taken together have the same
number of non-zero stencil entries as $A^{p}$. Therefore, we only have to
evaluate one additional stencil compared to the symmetric GS-Smoother. Note
that the _if_ statement inside the loop keeps the pseudo-algorithm concise and
in practice, it is avoided by moving it before the loop. Therefore, for two
sufficiently optimized implementations, we expect comparable computational
costs for the ILU and the symmetric GS-Smoother. In the following, we will
show that we can preserve the advantageous convergence rates from Section 3
with our surrogate smoother.

###  4.2 Numerical results

In the following, all numerical experiments are conducted on a mesh hierarchy
from levels 2 to 6. We use a V-cycle based multigrid solver with 3 pre- and
postsmoothing steps. A simple CG solver is used on the coarsest mesh. The
asymptotic convergence factor $\rho$ is determined by 20-steps of a power
iteration applied to the error propagation operator with a random initial
error.

Figure 6:  Convergence rates of the surrogate ILU for different degrees for
(V1) on the left and (V2) on the right.

To test the approximation quality of our ILU algorithm, we apply it to
different scalar coefficients $\kappa_{i}$, where

$\kappa_{i}(x,y,z)=1+10\left(x^{i}+y^{i}+z^{i}\right),\textmd{ with }0\leq
i\leq 3$

is a polynomial of degree $i$ on a trirectangular tetrahedron with unit height
and vertices $(0,0,0),(1,0,0),(0,1,0)$ and $(0,0,1)$. The asymptotic
convergence rates are given in Figure 6 for (V1) on the left and (V2) on the
right. The rates of the exact matrix version are plotted with dashed lines.
The operator $A$ is approximated by a surrogate with
$\textmd{dg}_{x}=\textmd{dg}_{y}=\textmd{dg}_{z}=3$ while the degree of the
surrogate approximation varies. For undistorted tetrahedra, the convergence
rates depend on the degree of the coefficient function $\kappa_{i}$. Using the
same degrees for the surrogate ILU and the coefficient function $\kappa_{i}$
allows us to recover the convergence rates of the original ILU algorithm. Both
versions of our algorithm perform similarly. Saving the boundary stencils in
(V1) only provides a small advantage in case of $\kappa_{3}$ compared to using
the surrogates everywhere and adjusting the boundary stencils in (V2), see
Figure 6.

Figure 7:  Convergence rates of the surrogate ILU for different degrees for
(V1) on the left and (V2) on the right.

In Figure 7, we use the constant $\kappa=1$ on a distorted tetrahedron of
height $h=0.1$. To approximate the Laplace operator for the residual, we use a
constant surrogate polynomial. On the left, we use isotropic degrees along
each coordinate axis for both variants of the surrogate ILU. The convergence
rates $\rho$ for both surrogate smoothers nearly coincide and saving the
boundary stencils yields only negligible improvements.

The previous results from Figure 3 suggest, that our stencil approximation has
to be accurate in the z-direction, while lower polynomial degrees should be
feasible inside the x-y-plane. Thus on the right of Figure 7, we set
$\textmd{dg}_{x}=\textmd{dg}_{y}=0$ and only vary $\textmd{dg}_{z}$. In this
case, during our iteration $\mathcal{O}(h^{-3})$ stencils do not have to be
calculated at all, only the $\mathcal{O}(h^{-1})$ stencils along the z-axis
have to be evaluated. The asymptotic convergence rate of the correct ILU are
not exactly attained even for large polynomial degrees, though the difference
is negligible for all practical purposes. To obtain a similar convergence rate
as for the isotropic case, a higher polynomial degree in z-direction has to be
used. Saving the boundary stencils leads here to a more robust version, though
again, both variants coincide if the polynomial degree in the z-direction is
large enough.

Figure 8:  Error in the discrete $L^{2}$-norm for the $L^{(\cdot)}_{d}$ and
$D^{(\cdot)}_{c}$ stencils of the distorted tetrahedron with height $h=0.1$
for levels $6$ and $7$. No values were skipped and no boundary stencils were
saved.

In Figure 8, the error in the discrete $L^{2}$-norm is plotted on different
multigrid levels for surrogate degrees from $0$ to $5$. As expected, for all
stencil directions, a higher polynomial degree leads to a smaller overall
error. For higher levels and hence finer meshes, this error becomes typically
larger for the same polynomial degrees. Recall that in Figure 4, we observed
that the non-asymptotic parts of the factorization are scaled in the spatial
direction by a factor of $\nicefrac{{1}}{{2}}$ from the coarser to the finer
level. Thus, the stencil functions which we aim to approximate obtain a larger
gradient and hence require higher polynomial degrees for the same
approximation error. An exception is the c-direction, whose absolute value
(Figure 4) decreases by a factor of $\nicefrac{{1}}{{2}}$ from the coarser
level to the finer level. This results in a smaller absolute error.

Figure 9:  Left: Test geometry with a distorted and blended tetrahedron.
Right: Convergence rates of the surrogate ILU on the given primitive for
different degrees.

To show that our surrogate ILU also works in scenarios with less artificial
geometries and coefficients, we apply it to a blended tetrahedron on the outer
boundary of a narrow spherical shell with outer radius 1 and inner radius 0.9
(Figure 9 left). The blended tetrahedron under consideration is depicted in
red while the unperturbed tetrahedron is indicated with black lines. The rest
of the shell which we do not use in the simulation is indicated in light-gray.
This benchmark combines both properties of the previous example: the varying
coefficient tensor due to the blending map and the strong anisotropy due to
the distorted base triangle.

For simplicity, we reassemble $A$ for each operator application instead of
relying on a surrogate approximation.

The convergence rates are depicted on the right of Figure 9, for surrogate
polynomials with isotropic degrees along the coordinate axes. Starting with
polynomial degree $3$, our ILU-Smoother is accurate enough to ensure
convergence. In this case, both variants of our ILU-Algorithm behave in the
same way. For polynomials of degrees larger than $7$, the rates for the
surrogate ILU and the matrix-based ILU coincide.

An ILU-factorization is not only used as a smoother for multigrid methods, but
can also be applied directly as a preconditioner for a different iterative
method. Candidates would be the Richardson iteration or Krylov methods. This
use-case is depicted in Figure 10: We apply the surrogate ILU as a
preconditioner within a CG method and increase its degree along the x-axis.

Figure 10:  Number of required iterations for a CG method directly
preconditioned by our ILU factorization to decrease the absolute residual
below $10^{-3}$ for the example from Fig. 9.

We present the required number of iterations for decreasing the absolute
residual below $10^{-3}$. The symmetric Gauss–Seidel method, not depicted in
the plot, needs $67$ iterations for achieving this. The ILU’s performance is
much better, starting at $41$ and $37$ iterations for low polynomial degrees
and resulting in $9$ to $7$ iterations for higher degrees. Similar to the
multigrid case, a surrogate polynomial degree of $7$ is roughly the point at
which our approximation quality stagnates and does not improve by much.

###  4.3 Performance analysis

| matrix-free | in-memory
---|---|---
line | memory traffic [Byte] | FLOP | memory traffic [Byte] | FLOP
5 | 0 | 0 | 0 | 0
6 | $15\cdot 8$ | 16 | $15\cdot 8$ | 16
11 | 0 | $7\cdot\textmd{deg}_{x}$ | $7\cdot 8$ | 0
12 | $6\cdot 8$ | 8 | $6\cdot 8$ | 8
$\sum$ | $21\cdot 8$ | $24+7\cdot\textmd{deg}_{x}$ | $28\cdot 8$ | 24
AI | $(\nicefrac{{1}}{{7}}+\nicefrac{{\textmd{deg}_{x}}}{{24}})$ FLOP/Byte | $\nicefrac{{3}}{{28}}$ FLOP/Byte
17 | 0 | $1\cdot\textmd{deg}_{x}$ | $1\cdot 8$ | 0
18 | $1\cdot 8$ | 1 | $1\cdot 8$ | 1
24 | 0 | $7\cdot\textmd{deg}_{x}$ | $7\cdot 8$ | 0
25 | $6\cdot 8$ | 8 | $6\cdot 8$ | 8
26 | $1\cdot 8$ | 1 | $1\cdot 8$ | 1
$\sum$ | $8\cdot 8$ | $10+8\cdot\textmd{deg}_{x}$ | $16\cdot 8$ | 10
AI | $(\nicefrac{{5}}{{32}}+\nicefrac{{\textmd{deg}_{x}}}{{8}})$ FLOP/Byte | $\nicefrac{{5}}{{64}}$ FLOP/Byte
Table 1: Theoretical cost analysis of memory traffic and FLOP in the interior
of our mesh.

Large scale low-order finite-element computations are typically memory-bound.
Therefore, to increase the performance of an algorithm we have to minimize
memory transfers. To judge the possible performance benefits of our surrogate
algorithm, we will therefore determine the arithmetic intensity of Algorithm 2
away from the boundary. The costs for the different lines of our
implementation are given in Table 1. The columns on the left relate to our
matrix-free implementation while the right column lists the theoretical costs
for a matrix-based implementation. We separately sum up the costs of the first
and second loop in Algorithm 2.

In this analysis, we assume that no blending and a constant coefficient is
used. Hence, the matrix stencil $A^{p}_{d}$ is constant, fits into the caches
and therefore Line 6 does not require memory traffic or FLOP. In Line 7, we
can assume that $x^{p+w}$ is already present in the cache. Therefore, only
$15$ vector entries of $x$ and $b$ have to be loaded, and $15$ multiplications
and one subtraction have to be applied. Since we assume to be in the interior,
Lines 7-10 do not contribute. In Line 11, evaluating each of the 7 stencil
entries takes $\textmd{deg}_{x}$ evaluations. In Line 12, we assume that
$w^{p}$ and $w^{p+w}$ are already cached and that 7 multiplications followed
by 1 addition are needed.

The in-memory implementation differs from the matrix-free implementation in
Lines 11 and 24, where instead of evaluating the stencil matrix-free, it has
to loaded from main memory. Altogether, $15\cdot 8$ additional bytes have to
be loaded from the main memory. Therefore, our matrix-free implementation has
the potential to decrease the necessary memory traffic by $34\%$.

Finally, we want to discuss the total memory requirements of the
factorization. Both for the matrix-free and the in-memory algorithm, we need
one additional intermediate vector $w$ which consists of $\textmd{\\#dof}\cdot
8$ bytes. Since the memory requirements for the polynomial coefficients can be
neglected, no additional memory for the matrix free implementation is
required. The in-memory algorithm on the other hand requires the storage of
the matrices $L$ and $D$. Since each interior row of $L$ consists of 7 non-
zeros, in total $\textmd{\\#dof}\cdot 64$ additional bytes have to be stored.
Therefore, the in-memory algorithm requires 9 times the memory of our matrix-
free implementation.

Figure 11:  Performance for the forward-substitution: Top: Runtime for
different surrogate polynomials. Bottom: Scalar double-precision FLOP/s for
different surrogate polynomials.

In Figure 11, we depict the experimentally determined runtime and the FLOP/s
averaged over the processes of the second variant of our forward substitution.
The code is executed on an Intel Xeon Gold 6136 with 12 processes distributed
over its two sockets. As an example, we used a grid with 48 macro-tetrahedra
on multigrid level 9 with 32 repetitions of the forward- and backward
substitutions, and as a measurement tool we rely on LIKWID [39, 44]. In the
figure, the performance of a simple matrix-based implementation, in which all
the stencils of the $L$, $D$ and $L^{T}$ matrices are stored continuously
inside a large memory buffer is depicted with a dashed line while the results
of an implementation which just relies on fixed asymptotic stencils is denoted
with a solid line. The latter provides a lower bound for the attainable
runtime of our algorithm, but has no practical relevance.

If we increase the polynomial degree in each direction uniformly, we arrive at
the  curve. With each increase in the degree, the FLOP count increases, and
brings a modest increase in runtime. For a degree larger or equal 3, storing
the stencils externally and loading them from memory is more efficient.

We showed that usually a high degree in the z-direction and lower degrees in
the x- and y-directions are completely sufficient to approximate the ILU-
Surrogate accurately enough. We use  and  to depict these scenarios, where the
degrees in x- and y-direction are fixed to 1 and 2, respectively, and only the
degree in z-direction is variable. In this case, the additional costs due to
the higher degrees in z-direction are negligible and are not visible in both
plots. The runtime for our surrogate based forward-substitution is
consistently lower than for the matrix-based version, even though we only get
a minor performance gain.

Overall, we have presented a memory efficient ILU-Smoother which can be
applied to much larger problems as if matrix-based implementations are used.
In the most general case, our method yields a modest runtime penalty for
larger polynomial degrees. However, due to the reduced memory traffic, we can
even achieve a small runtime gain with our surrogate approach for many
relevant cases.

##  5 Conclusion

In this paper, we have introduced an ILU algorithm on hybrid grid geometries
and investigated its performance within a multigrid solver. We replaced the
matrix-based algorithm by two matrix-free variants based on surrogate
polynomials approximating the stencils of the ILU matrix. To our knowledge,
this is the first matrix-free realization of a non-local operator based on an
algebraic factorization. Both matrix-free methods could attain the asymptotic
convergence rates of their matrix-based counterparts.

The ILU convergence rates were robust for individual distorted tetrahedra
which suggests a large performance gain for the matrix-free evaluation of the
Steklov–Poincaré operator on distorted hybrid grids compared to simpler
smoothing schemes as the GS algorithm. The surrogate ILU provides a huge
memory reduction with a small improvement in runtime compared to the usual
matrix-based realizations.

Furthermore, only the standard ILU-Smoother was investigated. The extension of
our algorithm to a matrix-free version of the modified ILU described in Remark
3.1 and an investigation of its performance are still of interest.

#### Code availability

The software used for obtaining the presented results are part of the open
source framework HyTeG [27] and publicly available at [1]. The results can be
reproduced by executing the Python scripts in apps/ILUSmoother/scripts.

##  Appendix A Smoother performance on single macro-tetrahedra

The following numerical experiments investigate the performance of our
smoother on single macro-tetrahedra.

Figure 12:  Tetrahedron shapes.

In Figure 12, we depict different tetrahedral shapes based on the ones in
[21], which are the “Spindle”, the “Cap”, the “Spade”, and a regular
tetrahedron. The red numbers are used to assign numbers to the vertices to
which a permutation will be applied. The tetrahedral coordinates are given in
Table SM1 in the supplementary materials.

Due to the symmetry of the shapes, we only obtain a small number of different
permutations. The asymptotic convergence rates for these are given in Table 2
for both our ILU and symmetric GS-Smoother. The rate determined by our
permutation heuristic is marked in bold. For all shapes, our ILU-Smoother
yields better rates than the GS-Smoother if the best orientation is used. This
gain is especially dramatic for the “Cap” shape, where
$\rho_{\textmd{GS}}\approx 0.5$, which means that 6-7 MG iterations are
necessary with GS when compared to a single MG iteration with an ILU-Smoother.
In addition, the optimal permutations for the GS and the ILU-Smoother always
coincide. For a fair comparison, we will therefore use in all our future
comparisons the same ordering for the GS as for the ILU-Smoother.

$\pi$ | $(1\,2\,3\,4)$ | $(1\,3\,2\,4)$ | $(1\,4\,2\,3)$ | $(1\,2\,3\,4)$ | $(1\,2\,4\,3)$ | $(1\,3\,4\,2)$ | $(2\,3\,4\,1)$
---|---|---|---|---|---|---|---
$\rho_{\textmd{GS}}$ | 0.77 | 0.54 | 0.78 | 0.52 | 0.53 | 0.52 | 0.51
$\rho_{\textmd{ILU}}$ | 0.65 | 0.39 | 0.35 | 0.010 | 0.43 | 0.43 | 0.0096
$\\#_{\textmd{GS}}$ | 53 | 23 | 56 | 22 | 22 | 22 | 21
$\\#_{\textmd{ILU}}$ | 33 | 15 | 14 | 3 | 17 | 17 | 3
| Spindle | Cap

$\pi$ | $(1\,2\,3\,4)$ | $(1\,2\,4\,3)$ | $(1\,3\,4\,2)$ | $(2\,1\,3\,4)$ | $(2\,1\,4\,3)$ | $(2\,3\,4\,1)$ | $(1\,2\,3\,4)$
---|---|---|---|---|---|---|---
$\rho_{\textmd{GS}}$ | 0.20 | 0.085 | 0.20 | 0.079 | 0.20 | 0.055 | 0.054
$\rho_{\textmd{ILU}}$ | 0.084 | 0.053 | 0.060 | 0.014 | 0.14 | 0.028 | 0.025
$\\#_{\textmd{GS}}$ | 9 | 6 | 9 | 6 | 9 | 5 | 5
$\\#_{\textmd{ILU}}$ | 6 | 5 | 5 | 4 | 8 | 4 | 4
| Spade | Regular

Table 2:  Asymptotic convergence rates $\rho$ and the number of iterations
$\\#$ to decrease the error by $10^{-6}$ for different permutations $\pi$. The
convergence rates due to our reordering algorithm are marked in bold.  Figure
13:  The multigrid algorithm on a distorted tetrahedron. Left: The asymptotic
convergence rates for the multigrid algorithm. Middle: The convergence rates
in a semilogarithmic plot. Right: Comparison of
$\log(\rho_{ILU})/\log(\rho_{SGS})$.

We depict the convergence rate of a trirectangular tetrahedron for different
heights of the top vertex in Figure 13. As we see in the leftmost plot,
decreasing the height degenerates the asymptotic convergence rate $\rho$ for
multigrid with symmetric GS smoothing. It is a well-known fact that the ILU-
Smoother remains robust with respect to these deformations. The logarithmic
plot in the middle even suggests, that it becomes an exact solver. In the
right we have plotted $\log(\rho_{ILU})/\log(\rho_{SGS})$, which is a measure
how many additional multigrid iterations the SGS based multigrid need to
achieve the same convergence rate as the ILU based multigrid. The performance
gain due to an ILU becomes arbitrary large for small heights $h$.

Both Table 2 and Figure 13 suggest, that the ILU on the reoriented tetrahedra
mesh provides a robust smoother inside the macro-tetrahedrons. This transfers
directly to the Steklov–Poincaré operator in which we use our ILU as an inner
solver component.

##  Appendix SM1 Supplementary Material: Hybrid-Smoother

In this second application, we extend the subgrid ILU-Smoother to a smoother
on the global grid. On the interfaces between the structured grids, we apply a
simple Gauss–Seidel-Smoother (GS-Smoother). Thereby, we obtain a block
smoother that we examine on a hybrid structured tetrahedral grid for a finite-
element discretization of piecewise continuous finite-element functions.

Given a set of preconditioner matrices
$C_{p}\in\mathbb{R}^{|\mathcal{I}_{p}|\times|\mathcal{I}_{p}|}$ for the lower
dimensional primitives
$p\in\mathcal{V}_{H}\cup\mathcal{E}_{H}\cup\mathcal{F}_{H}$ and a set of
symmetric preconditioner matrices
$C_{t,sym}=C^{T}_{t,sym}\in\mathbb{R}^{|\mathcal{I}_{t}|\times|\mathcal{I}_{t}|}$
for each tetrahedron $t\in\mathcal{T}_{H}$, we define a symmetric hybrid
smoother in Algorithm 3.

1:$\mathbf{u}\leftarrow\mathbf{u}+\sum_{v\in\mathcal{V}_{H}}R_{v}^{T}C^{-1}_{v}R_{v}\,(\mathbf{f}-A\mathbf{u})$.

2:$\mathbf{u}\leftarrow\mathbf{u}+\sum_{e\in\mathcal{E}_{H}}R_{e}^{T}C^{-1}_{e}R_{e}\;(\mathbf{f}-A\mathbf{u})$.

3:$\mathbf{u}\leftarrow\mathbf{u}+\sum_{f\in\mathcal{F}_{H}}R_{f}^{T}C^{-1}_{f}R_{f}(\mathbf{f}-A\mathbf{u})$.

4:$\mathbf{u}\leftarrow\mathbf{u}+\sum_{t\in\mathcal{T}_{H}}R_{t}^{T}C^{-1}_{t,sym}R_{t}(\mathbf{f}-A\mathbf{u})$.

5:$\mathbf{u}\leftarrow\mathbf{u}+\sum_{f\in\mathcal{F}_{H}}R_{f}^{T}C^{-T}_{f}R_{f}(\mathbf{f}-A\mathbf{u})$.

6:$\mathbf{u}\leftarrow\mathbf{u}+\sum_{e\in\mathcal{E}_{H}}R_{e}^{T}C^{-T}_{e}R_{e}\;(\mathbf{f}-A\mathbf{u})$.

7:$\mathbf{u}\leftarrow\mathbf{u}+\sum_{v\in\mathcal{V}_{H}}R_{v}^{T}C^{-T}_{v}R_{v}\,(\mathbf{f}-A\mathbf{u})$.

Algorithm 3 Symmetric hybrid smoother.

This smoother sequentially updates the degrees-of-freedom (DoF) on the macro-
vertices, macro-edges, macro-faces and macro-tetrahedra. In order to retain
symmetry, we traverse the macro-hierarchy in the reverse direction by applying
smoothing steps to the macro-face, macro-edge and finally to the macro-vertex
DoF. Note that procedures like that can be performed efficiently for
distributed memory parallelizations since all DoF sharing a primitive type can
be smoothed in parallel [27]. Due to the symmetrization it is possible to use
our multigrid algorithm as a preconditioner for a conjugate gradient (CG) or
minimial residual (MINRES) method.

We define a hybrid symmetric Gauss–Seidel-Smoother (SGS-Smoother) as a
reference example. It is given by

$C_{v}=L_{A,v},\,C_{e}=L_{A,e},\,C_{f}=L_{A,f},\textmd{ and
}C_{t,sym}=(L_{A,t})(D_{A,t})^{-1}(L_{A,t})^{T},$

for $v\in\mathcal{V}_{H}$, $e\in\mathcal{E}_{H}$, $f\in\mathcal{F}_{H}$ and
$t\in\mathcal{T}_{H}$. For simplicity, we will refer to it by Gauss–Seidel
even though this is technically not correct since some of the DoF on the
interfaces $\mathcal{V}_{H}\cup\mathcal{E}_{H}\cup\mathcal{F}_{H}$ are handled
additively instead of multiplicatively. Numerical experiments in [27] have
shown that in practice this has no impact on the convergence rates.

The hybrid ILU-Smoother is defined similarly by

$C_{v}=L_{A,v},\,C_{e}=L_{A,e},\,C_{f}=L_{A,f},\quad\textmd{ and }\quad
C_{t,sym}=L_{t}D_{t}L_{t}^{T}.$

In a nutshell, the smoother consists of Gauss–Seidel steps on the macro-
vertices, macro-edges, and macro-faces. The ILU itself is just applied inside
the macro-tetrahedra.

######  Remark SM1.1.

We stress that the submatrices $R_{v}AR^{T}_{v}$, $R_{e}AR^{T}_{e}$, and
$R_{f}AR^{T}_{f}$ for $v\in\mathcal{V}_{H},e\in\mathcal{E}_{H},$ and
$f\in\mathcal{F}_{H}$ on the lower dimensional primitives are well conditioned
matrices. Hence, even a simple iterative scheme is sufficient for reducing the
error significantly. Therefore, using an ILU also on the lower primitives is
not necessary and only complicates the implementation.

###  SM1.1 Numerical experiments

We investigate the performance of our hybrid ILU-Smoother on a hybrid mesh
consisting of several tetrahedra. For this, we will introduce a suitable
benchmark problem. Firstly, we will directly apply the multigrid algorithm as
solver and in a second scenario we will use the multigrid as a preconditioner
within a CG method.

#### Benchmark

Figure 1 depicts our benchmark scenario. A unit cube with Dirichlet boundary
conditions on the top and bottom and Neumann boundary conditions on the
remaining sides is divided into two parts where $h_{\textmd{lower}}$ is the
height of the lower volume. In both volumes, we assume a piecewise constant
material parameter $\kappa$ with a possible jump at the interface. We want to
compare the convergence rates for different heights of the lower cube. This is
similar to a scenario in geophysics in which the viscosity jumps from a narrow
upper layer in the Earth’s mantle that can barely be resolved to a wider lower
layer [16]. We will use this scenario to assess the performance of our hybrid
ILU-Smoother.

Fig. SM 1:  Benchmark scenario: The size of the lower half of the cube is
decreased to approximate a tiny layer in the Earth’s mantle.

Geometry | x | y | z
---|---|---|---
Spindle | 0.0 | 0.0 | 0.5
| 0.0 | 0.0 | -0.5
| 0.5 | 1.0 | 0.0
| -0.5 | 1.0 | 0.0
Cap | 0.0 | 0.0 | 0.0
| 1.0 | 0.0 | 0.0
| 0.5 | 0.866 | 0.0
| 0.5 | 0.288 | 0.093
Spade | 0.0 | 0.0 | 0.0
| 1.0 | -0.666 | 0.0
| 1.0 | 0.666 | 0.0
| 1.0 | 0.0 | 0.443

Table SM 1:  Tetrahedral coordinates for the Spindle, Cap and Spade.

#### Multigrid solver

The results for the asymptotic convergence rates of our multigrid algorithm
are depicted in the left of Figure 2, for a fixed $\kappa_{lower}=1$ and 3
different choices of $\kappa_{upper}$. As expected, since the material jump
happens on an interface that can be resolved by all coarse grids, it does not
impact the SGS-Smoother (see [23, Sec. 10.3]) and the same is true for our
ILU-Smoother. For jumps inside the elements, special interpolation operators
would have to be introduced [2]. Our ILU-Smoother consistently performs better
than the SGS-Smoother for the multigrid algorithm. For small tetrahedral
heights, the performance of both smoothers decays and the convergence rate
approaches one. The slope of the performance degradation is similar for both
smoothers. Note that on the interfaces, we use the same Gauss–Seidel smoothing
strategy for both smoothers which suggests that the performance loss is
related to them and better methods are necessary to handle the interfaces.

Fig. SM 2:  Different solvers applied to our benchmark scenario (Fig. 1).
Left: Asymptotic convergence rates with a multigrid V-Cycle on grid levels 2
to 6, and 3 pre- and postsmoothing steps. Right: The V-Cycle used as a
preconditioner for a PCG solver with iteration numbers in the upper plot and
the ratio of the SGS and ILU iteration numbers in the lower plot.

#### Multigrid preconditioner

In the right plot of Figure 2, we use our multigrid algorithm only as a
preconditioner for a CG solver and plot the iteration count for getting the
absolute unpreconditioned residual below $10^{-5}$. The hybrid ILU-Smoother is
an improvement compared to the SGS-Smoother, but again the hybrid solver
cannot preserve the robustness with respect to degenerated tetrahedra. The
ratio between the number of iterations between a GS and an ILU based multigrid
preconditioner is depicted in the lower right. For a large range of heights it
stays around 1.2 before it degenerates to 1 for extremely small heights.

###  SM1.2 Conclusion

We extended the smoother to a hybrid ILU-Smoother to directly use it as a
smoother on a more complex domain with several tetrahedrons joined via an
interface. In this case, our algorithm outperforms the SGS-Smoother while
having comparable computational costs.

Although, a hybrid block smoother is perfectly suited for matrix-free
approaches on hybrid meshes, we observed that it lacks robustness with respect
to poorly structured macro-meshes. Future work might therefore include
applying more sophisticated techniques on the interfaces to further improve
the performance.

##  Appendix SM2 Simulation data

Table 1 provides the detailed coordinates of the reference tetrahedra for
result reproductions and comparisons.

## References

  * [1] ILU implementation - source code.  https://doi.org/10.5281/zenodo.7199022.
  * [2] R. E. Alcouffe, A. Brandt, J. E. Dendy, Jr, and J. W. Painter, The multi-grid method for the diffusion equation with strongly discontinuous coefficients, SIAM Journal on Scientific and Statistical Computing, 2 (1981), pp. 430–454.
  * [3] H. Anzt, E. Chow, and J. Dongarra, ParILUT—a new parallel threshold ILU factorization, SIAM Journal on Scientific Computing, 40 (2018), pp. C503–C519.
  * [4] H. Anzt, T. Ribizel, G. Flegar, E. Chow, and J. Dongarra, ParILUT-a parallel threshold ILU for GPUs, in 2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS), IEEE, 2019, pp. 231–241.
  * [5] O. Axelsson, Incomplete block matrix factorization preconditioning methods. the ultimate answer?, Journal of Computational and Applied Mathematics, 12 (1985), pp. 3–18.
  * [6] S. Bauer, H.-P. Bunge, D. Drzisga, S. Ghelichkhan, M. Huber, N. Kohl, M. Mohr, U. Rüde, D. Thönnes, and B. Wohlmuth, TerraNeo—mantle convection beyond a trillion degrees of freedom, in Software for Exascale Computing-SPPEXA 2016-2019, Springer, 2020, pp. 569–610.
  * [7] S. Bauer, D. Drzisga, M. Mohr, U. Rüde, C. Waluga, and B. Wohlmuth, A stencil scaling approach for accelerating matrix-free finite element implementations, SIAM Journal on Scientific Computing, 40 (2018), pp. C748–C778.
  * [8] S. Bauer, M. Huber, S. Ghelichkhan, M. Mohr, U. Rüde, and B. Wohlmuth, Large-scale simulation of mantle convection based on a new matrix-free approach, Journal of Computational Science, 31 (2019), pp. 60–76, https://doi.org/10.1016/j.jocs.2018.12.006, https://doi.org/10.1016/j.jocs.2018.12.006.
  * [9] S. Bauer, M. Huber, M. Mohr, U. Rüde, and B. Wohlmuth, A new matrix-free approach for large-scale geodynamic simulations and its performance, in International Conference on Computational Science, Springer, 2018, pp. 17–30.
  * [10] S. Bauer, M. Mohr, U. Rüde, J. Weismüller, M. Wittmann, and B. Wohlmuth, A two-scale approach for efficient on-the-fly operator assembly in massively parallel high performance multigrid codes, Applied Numerical Mathematics, 122 (2017), pp. 14 – 38, https://doi.org/https://doi.org/10.1016/j.apnum.2017.07.006, http://www.sciencedirect.com/science/article/pii/S0168927417301642.
  * [11] B. Bergen, Hierarchical Hybrid Grids: Data Structures and Core Algorithms for Efficient Finite Element Simulations on Supercomputers, SCS Publishing House, Erlangen, 2005.
  * [12] B. Bergen and F. Hülsemann, Hierarchical hybrid grids: data structures and core algorithms for multigrid, Numer. Lin. Alg. Appl., 11 (2004), pp. 279–291.
  * [13] B. Bergen, G. Wellein, F. Hülsemann, and U. Rüde, Hierarchical hybrid grids: Achieving TERAFLOP performance on large scale finite element simulations, International Journal of Parallel, Emergent and Distributed Systems, 22 (2007), pp. 311–329.
  * [14] J. Bey, Tetrahedral grid refinement, Computing, 55 (1995), pp. 355–378.
  * [15] T. F. Chan and T. P. Mathew, Domain decomposition algorithms, Acta numerica, 3 (1994), pp. 61–143.
  * [16] G. F. Davies and M. A. Richards, Mantle convection, The Journal of Geology, 100 (1992), pp. 151–206.
  * [17] D. Drzisga, B. Keith, and B. Wohlmuth, The surrogate matrix methodology: A priori error estimation, SIAM Journal on Scientific Computing, 41 (2019), pp. A3806–A3838, https://doi.org/10.1137/18M1226580.
  * [18] D. Drzisga, B. Keith, and B. Wohlmuth, The surrogate matrix methodology: Low-cost assembly for isogeometric analysis, Computer Methods in Applied Mechanics and Engineering, 361 (2020), p. 112776, https://doi.org/10.1016/j.cma.2019.112776.
  * [19] D. Drzisga, U. Rüde, and B. Wohlmuth, Stencil scaling for vector-valued PDEs on hybrid grids with applications to generalized Newtonian fluids, SIAM Journal on Scientific Computing, 42 (2020), pp. B1429–B1461.
  * [20] F. J. Gaspar, J. L. Gracia, and F. J. Lisbona, Fourier analysis for multigrid methods on triangular grids, SIAM Journal on Scientific Computing, 31 (2009), pp. 2081–2102.
  * [21] B. Gmeiner, T. Gradl, F. Gaspar, and U. Rüde, Optimization of the multigrid-convergence rate on semi-structured meshes by local Fourier analysis, Computers & Mathematics with Applications, 65 (2013), pp. 694–711.
  * [22] L. Grigori and S. Moufawad, Communication avoiding ILU0 preconditioner, SIAM Journal on Scientific Computing, 37 (2015), pp. C217–C246.
  * [23] W. Hackbusch, Multi-grid methods and applications, vol. 4, Springer Science & Business Media, 2013.
  * [24] S. Kang, L. C. Ngo, H. Choi, W. Chung, Y.-H. Yoo, and J. Y. Yoo, Performance comparison of parallel ILU preconditioners for the incompressible Navier-Stokes equations, Journal of Mechanical Science and Technology, 34 (2020), pp. 1175–1184.
  * [25] G. Karypis and V. Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM Journal on Scientific Computing, 20 (1998), pp. 359–392, https://doi.org/10.1137/S1064827595287997, https://doi.org/10.1137/S1064827595287997, https://arxiv.org/abs/https://doi.org/10.1137/S1064827595287997.
  * [26] R. Kettler and P. Wesseling, Aspects of multigrid methods for problems in three dimensions, Applied mathematics and computation, 19 (1986), pp. 159–168.
  * [27] N. Kohl, D. Thönnes, D. Drzisga, D. Bartuschat, and U. Rüde, The HyTeG finite-element software framework for scalable multigrid solvers, International Journal of Parallel, Emergent and Distributed Systems, 34 (2019), pp. 477–496.
  * [28] V. G. Korneev and U. Langer, Domain Decomposition Methods and Preconditioning, John Wiley & Sons, Ltd, 2004, ch. 22, https://doi.org/https://doi.org/10.1002/0470091355.ecm019, https://onlinelibrary.wiley.com/doi/abs/10.1002/0470091355.ecm019, https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/0470091355.ecm019.
  * [29] M. Kronbichler, K. Ljungkvist, M. Allalen, M. Ohlerich, I. Pasichnyk, and W. A. Wall, Performance optimization of matrix-free finite-element algorithms within deal. II, (2017).
  * [30] D. Lukarski, H. Anzt, S. Tomov, and J. Dongarra, Hybrid multi-elimination ILU preconditioners on GPUs, in 2014 IEEE International Parallel Distributed Processing Symposium Workshops, 2014, pp. 7–16, https://doi.org/10.1109/IPDPSW.2014.7.
  * [31] T. Malas and L. Gürel, Incomplete LU preconditioning with the multilevel fast multipole algorithm for electromagnetic scattering, SIAM Journal on Scientific Computing, 29 (2007), pp. 1476–1494.
  * [32] T. Mathew, Domain decomposition methods for the numerical solution of partial differential equations, vol. 61, Springer Science & Business Media, 2008.
  * [33] M. Mayr, L. Berger-Vergiat, P. Ohm, and R. S. Tuminaro, Non-invasive multigrid for semi-structured grids, arXiv preprint arXiv:2103.11962, (2021).
  * [34] J. A. Meijerink and H. A. Van Der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric $m$-matrix, Mathematics of computation, 31 (1977), pp. 148–162.
  * [35] K.-D. Oertel and K. Stüben, Multigrid with ILU-smoothing: Systematic tests and improvements, in Robust Multi-Grid Methods, Springer, 1989, pp. 188–199.
  * [36] F. Pellegrini, Scotch and PT-Scotch Graph Partitioning Software: An Overview, in Combinatorial Scientific Computing, O. S. Uwe Naumann, ed., Chapman and Hall/CRC, 2012, pp. 373–406, https://doi.org/10.1201/b11644-15, https://hal.inria.fr/hal-00770422.
  * [37] M. Pinto, C. Rodrigo, F. Gaspar, and C. Oosterlee, On the robustness of ILU smoothers on triangular grids, Applied Numerical Mathematics, 106 (2016), pp. 37–52.
  * [38] A. Quarteroni and A. Valli, Domain decomposition methods for partial differential equations, no. BOOK, Oxford University Press, 1999.
  * [39] T. Roehl, J. Treibig, G. Hager, and G. Wellein, Overhead analysis of performance counter measurements, in 43rd International Conference on Parallel Processing Workshops (ICCPW), Sept 2014, pp. 176–185, https://doi.org/10.1109/ICPPW.2014.34.
  * [40] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
  * [41] R. Tielen, M. Möller, D. Göddeke, and C. Vuik, p-multigrid methods and their comparison to h-multigrid methods within isogeometric analysis, Computer Methods in Applied Mechanics and Engineering, 372 (2020), p. 113347.
  * [42] R. Tielen, M. Möller, and K. Vuik, A direct projection to low-order level for p-multigrid methods in isogeometric analysis, in Numerical Mathematics and Advanced Applications ENUMATH 2019, Springer, 2021, pp. 1001–1009.
  * [43] A. Toselli and O. Widlund, Domain decomposition methods-algorithms and theory, vol. 34, Springer Science & Business Media, 2004.
  * [44] J. Treibig, G. Hager, and G. Wellein, LIKWID: A lightweight performance-oriented tool suite for x86 multicore environments, in Proceedings of PSTI2010, the First International Workshop on Parallel Software Tools and Tool Infrastructures, San Diego CA, 2010.
  * [45] N. Umetani, S. P. MacLachlan, and C. W. Oosterlee, A multigrid-based shifted Laplacian preconditioner for a fourth-order Helmholtz discretization, Numerical Linear Algebra with Applications, 16 (2009), pp. 603–626.
  * [46] A. Vargas, T. M. Stitt, K. Weiss, V. Z. Tomov, J.-S. Camier, T. Kolev, and R. N. Rieben, Matrix-free approaches for GPU acceleration of a high-order finite element hydrodynamics application using MFEM, Umpire, and RAJA, arXiv preprint arXiv:2112.07075, (2021).
  * [47] P. Wesseling, Introduction to multigrid methods, tech. report, 1995\.