---
abstract: 'Local Fourier analysis is a commonly used tool to assess the quality and aid in the construction of geometric multigrid methods for translationally invariant operators. In this paper we automate the process of local Fourier analysis and present a framework that can be applied to arbitrary, including non-orthogonal, repetitive structures. To this end we introduce the notion of crystal structures and a suitable definition of corresponding wave functions, which allow for a natural representation of almost all translationally invariant operators that are encountered in applications, e.g., discretizations of systems of PDEs, tight-binding Hamiltonians of crystalline structures, colored domain decomposition approaches and last but not least two- or multigrid hierarchies. Based on this definition we are able to automate the process of local Fourier analysis both with respect to spatial manipulations of operators as well as the Fourier analysis back-end. This automation most notably simplifies the user input by removing the necessity for compatible representations of the involved operators. Each individual operator and its corresponding structure can be provided in any representation chosen by the user.'
author:
- 'Karsten Kahl[^1]'
- Nils Kintscher
bibliography:
- 'Manuscript.bib'
title: 'Automated Local Fourier Analysis (aLFA)[^2]'
---

Introduction
============

Local Fourier analysis (LFA) is a powerful tool used in the construction and analysis of multigrid methods introduced in [@Bran1977]. The fundamental idea of LFA is to leverage the connection between position space and frequency space via the Fourier transform. That is, in case the involved operators can be described by stencils in position space, meaning that they are translationally invariant, their Fourier transform yields so-called symbols, which can be handled much more easily. In the context of multigrid methods LFA can be used to obtain precise approximations of the asymptotic convergence rate by assessing the spectral radius of the corresponding error propagation operator [@Saad]. This approximation of the convergence rate is (asymptotically) still valid if the translational invariance is slightly violated, as for example when lexicographic Gauss-Seidel type smoothers are used, or in the case of certain non-periodic boundary conditions [@Bran1994; @MR3866073; @Stevenson1990]. In other cases a similar convergence rate can usually be obtained by additional processing [@KahlKint2018; @TrotOostSchu2001]. Due to these facts, LFA is one of the main tools in the quantitative analysis of two- and multigrid methods. An introduction to LFA including several examples can be found in [@TrotOostSchu2001; @WienJopp2004]. Multigrid methods have first been considered for the solution of the linear systems of equations originating in the discretization of (elliptic) partial differential equations (PDEs) [@TrotOostSchu2001]. Due to the fact that the simplest tiling of space is rectangular and discretizations are particularly simple to carry out on such tilings, the usual multigrid components and the LFA have originally been designed and tailored for such discretizations (cf. [@TrotOostSchu2001; @WienJopp2004]). Several other geometries, including systems of PDEs, have been considered in the past as well. LFA has been carried out for operators defined on triangular tilings in [@GaspGracLisb2009] and on hexagonal tilings in [@ZhouFult]. Further, it has been applied to edge-based quadrilateral discretizations [@BoonLentVand2008], regular Voronoi meshes associated with acute triangular grids [@RodrSaliGaspLisb2014], edge-based discretizations on triangular grids [@rodrigo_sanz_gaspar_lisbona_2015] and jumping coefficients on rectangular grids [@RitExtending2017; @RittichBolten2018]. These papers do a complete two-grid analysis, and in some cases even a three-grid analysis, which was first introduced in [@WienOost2001].

While the concept of LFA is well understood, its application quickly becomes complex and involved the more frequencies get intermixed, e.g., by block smoothers, in a three-grid analysis, or in higher dimensional problems ($n > 2$). Thus, there exists software that automates the application of the LFA [@lfalab_HR; @WienJopp2004]. In contrast to the software described in [@WienJopp2004], which basically consists of a collection of templates corresponding to certain smoother and restriction/prolongation strategies for specific problems, the software [@lfalab_HR], freely available on GitHub,[^3] can be used to analyze arbitrary translationally invariant operators on rectangular grids. This software has for example been used to analyze colored block Jacobi methods in combination with aggressive coarsening applied to PDEs with jumping coefficients in [@RitExtending2017; @RittichBolten2018]. As the number of frequencies which get intermixed increases with the block-size of the smoother and the growing coarsening rate, a manual analysis of this problem would be laborious. In contrast to the approach developed in this paper, the bulk of the analysis in the previously mentioned references is mainly carried out in frequency space. LFA with an emphasize on position space is rare in the literature. In [@Huckle2008] under the name *compact Fourier analysis*, a position space oriented LFA is introduced where block Toeplitz matrices are used in order to capture the different classes of unknowns. This work has some aspects in common with the approach to LFA we develop in this paper, but lacks generality as it is limited to simple systems on rectangular grids.

The LFA presented in this paper unifies the position space oriented approaches and allows the treatment of operators on arbitrary repetitive structures. To do so we introduce a mathematical framework for the analysis of translationally invariant operators which alter value distributions on lattices and crystals [@Ashcroft; @Schrijver:1986:TLI:17634]. These structures can be based on arbitrary sets of primitive vectors, including non-orthogonal ones, i.e., non-rectangular structures. These crystal structures, which naturally occur in the tight-binding descriptions of solid-state physics [@Ashcroft] or discretizations of systems of PDEs, enable the convenient and concise description of the resulting operators and allow for the automatic generation of their representations when enlarging their translational invariance, e.g., coarsening in the context of multigrid methods. Furthermore, they are very helpful in the representation of overlapping and non-overlapping block smoothers. Our framework is developed to such an extent that the only task required of the user is to provide a description of the occurring operators with respect to (potentially non-matching) descriptions of the underlying repetitive structures, i.e., each operator can be supplied in the simplest or most convenient representation. The remainder of the analysis can then be carried out automatically. In contrast to previously developed LFA, this is achieved by explicitly including a connection of the operator to its underlying structure. This allows us on one hand to automate the transformation of operators in position space, e.g., by finding a least common lattice of translational invariance of two operators and to rewrite their representations accordingly. On the other hand, this focus on structure yields a natural representation and discretization of the dual space that enables the automation of the frequency space part of the analysis as well. All these tasks can be carried out using basic principles and normal forms of integer linear algebra [@Schrijver:1986:TLI:17634]. In combination these developments alleviate the use of LFA by removing any manual calculations. That is, neither a mixing analysis nor transformations of operators to common (and rectangular) translational invariances have to be carried out by hand. While our automated LFA does not necessarily enlarge the set of methods that are analyzable by conventional LFA, it enables the reliable and easy-to-use analysis of complex methods on complicated structures (e.g., overlapping block smoothers and discretizations of systems of PDEs). An open-source Python implementation of the automated LFA framework [@aLFA_NK] is freely available on GitLab.[^4]

The automation presented in this paper does have some limitation in terms of the smoothers that can be analyzed. Any sequential, i.e., lexicographic, smoother with overlapping update regions changes values in the overlap multiple times in one application. This cannot be easily translated to the structures introduced in this paper. While such smoothers have been analyzed in frequency space before (cf. [@MacLOost2011; @Molenaar1991; @Sivaloganathan:1991:ULM]), this particular treatment of sequential overlap is momentarily not covered in our framework. Note, that the mere presence of overlap is not the problem here. By introducing a coloring, such that the complete sweep can be split into a sequence of updates where each one of them only changes values at most once, automated LFA can be applied as we are going to demonstrate in this paper. Due to the fact that a coloring in overlapping approaches also favors parallelism over their sequential counterparts, we feel that this limitation is relatively minor when targeting actual applications.

The paper is organized as follows. In \[sec:lattices\_crystals\] we introduce basic notation to describe the underlying structures of the operators we want to analyze: lattices and crystals. In the context of LFA we are interested in translationally invariant operators which alter value distributions on crystals. These kind of operators and their properties are specified in \[sec:operators\]. After that, \[sec:example\] illustrates the introduced notation by means of the discretized Laplacian and the red-black Gauss-Seidel method, which is the simplest method where a difference to the conventional LFA becomes apparent. In here, we manually rewrite the discretized Laplacian with respect to the translational invariance of the red-black splitting, such that the method can be analyzed via the results of \[sec:operators\]. contains the general theoretical results which are used to automate the complete procedure, removing the need to manually modify operators. In here, several results of integer linear algebra are used, which we review in \[sec:sublat\_quotient\]. Using these results we obtain the algorithms given in \[sec:algorithms\] which, in combination with the arithmetic of multiplication operators given in \[sec:rules\_of\_comp\], realize the automated LFA. Finally, \[sec:Application\] contains selected examples to demonstrate the merits of the developed approach. First, we analyze a $4$-color overlapping block Gauss-Seidel smoother for the tight-binding Hamiltonian of the carbon allotrope graphene and give some two-grid convergence results. Second, we reproduce the two-level analysis for the curl-curl problem found in [@BoonLentVand2008] as a further illustration of how our approach is applied to complex error propagators and to double-check its results.

Lattices and Crystals {#sec:lattices_crystals}
=====================

In order to be able to automate the process of LFA within a unified framework we first have to review some basic definitions of integer linear algebra and crystallography [@Ashcroft; @Schrijver:1986:TLI:17634]. An (ideal) crystal is an infinite repetition, defined by a lattice, of a structure element.

Let ${\mathcal{a}}_1,{\mathcal{a}}_2,\ldots,{\mathcal{a}}_n \in {
    \mathbb{R}}^n$ be linearly independent. An $n$-dimensional lattice ${
    \mathbb{L}}$ is the set of points $$\begin{aligned}
            {
    \mathbb{L}}={\{x=\sum_{\ell=1}^n j_\ell {\mathcal{a}}_\ell \in {
    \mathbb{R}}^n \, : \, j_1, j_2,\ldots, j_n \in {
    \mathbb{Z}}\}}.
        \end{aligned}$$ The vectors ${\mathcal{a}}_1,{\mathcal{a}}_2,\ldots,{\mathcal{a}}_n$ are known as the *primitive vectors* or *lattice basis*. Using matrix notation, i.e., ${\mathcal{A}}:=\mat[c]{{\mathcal{a}}_1&{\mathcal{a}}_2&\ldots&{\mathcal{a}}_n}$, we can abbreviate the notation by $${
    \mathbb{L}}({\mathcal{A}}) :=  {\mathcal{A}}{
    \mathbb{Z}}^n  = {
    \mathbb{L}}.$$

Without loss of generality we restrict the definition of the second component of a crystal, the structure element, to primitive cells of the lattice.

\[def:primitive\_cell\] A *primitive cell* ${\Xi}={\Xi}({\mathcal{A}}) \subset {
    \mathbb{R}}^n$ of a lattice ${
    \mathbb{L}}({\mathcal{A}})$ is a (connected) volume of space that, if translated by all vectors of ${
    \mathbb{L}}$, fills up ${
    \mathbb{R}}^n$ completely without any overlap, i.e., $$\dot\cup_{x\in{
    \mathbb{L}}} {\{x + \xi  \, : \, \xi \in {\Xi}\}} ={
    \mathbb{R}}^n.$$ A common choice of primitive cells is given by $${{
    \mathcal{P}}}({\mathcal{A}}):={\mathcal{A}}[0,1)^n={\{y\in{
    \mathbb{R}}^n \, : \,  y = \sum_{\ell=1}^n \alpha_\ell {\mathcal{a}}_\ell,\, 0 \leq \alpha_\ell < 1\}},$$ i.e., parallelotopes spanned by the primitive vectors of ${
    \mathbb{L}}({\mathcal{A}})$.

The structure element of a crystal is now defined to consist of points ${\mathfrak{s}}_{1},\ldots,{\mathfrak{s}}_{m}$ that are contained in a particular primitive cell ${\Xi}$.

Let ${
    \mathbb{L}}({\mathcal{A}})$ be a lattice and ${\mathfrak{s}}\in {\Xi}({\mathcal{A}})^m$, $m \in \mathbb{N}$ be the *structure element*. A *crystal* is defined as the set of tuples $$\begin{aligned}
            {
    \mathbb{L}}^{{\mathfrak{s}}}({\mathcal{A}}):= {\{(x+{\mathfrak{s}}_1,x+{\mathfrak{s}}_2,\ldots,x+{\mathfrak{s}}_m) \, : \, x\in{
    \mathbb{L}}({\mathcal{A}}),\  {\mathfrak{s}}=({\mathfrak{s}}_1,\ldots,{\mathfrak{s}}_m)\}}.
        \end{aligned}$$ The elements of ${
    \mathbb{L}}^{{\mathfrak{s}}}({\mathcal{A}})$ are collectively written as $x+{\mathfrak{s}}=(x+{\mathfrak{s}}_1,x+{\mathfrak{s}}_2,\ldots,x+{\mathfrak{s}}_m)$.

We define a crystal and the associated structure element to be a tuple instead of a set as we want to study value distributions on crystals and particular operators which manipulate them. For this purpose, the order of a structure element is of importance.

To give an idea of the various occurrences of crystal structures in numerical applications we illustrate typical examples in \[fig:crystal\_examples\] together with their crystal representation. There are three main sources of repetitive structures that are well suited for crystal representations. First, the discretization of systems of PDEs on lattices lead to crystal structures, where the different species of unknowns typically constitute the structure element. Second, tight-binding Hamiltonian formulations in solid state physics for crystalline materials naturally imply a crystal representation based on the atomic structure. Last, colored domain decomposition approaches (e.g., red-black Gauss-Seidel) can be easily represented using crystals, where the smallest structure element typically consists of the union of one domain of each color.

Sublattices and quotient spaces {#sec:sublat_quotient}
-------------------------------

There are infinitely many representations of a crystal ${
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$. On one hand, the representation of any lattice in $n>1$ dimensions is non-unique, i.e., there exist different sets of primitive vectors that yield the same lattice structure. On the other hand, different representations of the crystal can be obtained by shifting the structure element or manipulating the underlying lattice structure, e.g., by using integer linear combinations of the primitive vectors and adjusting the structure element accordingly.

In order to cope with this lack of uniqueness of the representation of crystals, illustrated in \[fig:crystal\_example\], we introduce basic results from integer linear algebra [@Schrijver:1986:TLI:17634], which resolve the relationship of lattice structures.[^5] An important tool in this is the notion of a sublattice.

Let ${
    \mathbb{L}}({\mathcal{A}})$ and ${
    \mathbb{L}}({\mathcal{C}})$ be two lattices, if ${
    \mathbb{L}}({\mathcal{A}}) \supset {
    \mathbb{L}}({\mathcal{C}})$ then ${
    \mathbb{L}}({\mathcal{C}})$ is called *sublattice* of ${
    \mathbb{L}}({\mathcal{A}})$.

\[lem:sublattice\] A lattice ${
    \mathbb{L}}({\mathcal{C}})$ is a sublattice of ${
    \mathbb{L}}({\mathcal{A}})$ if and only if ${\mathcal{A}}^{-1} {\mathcal{C}}\in {
    \mathbb{Z}}^{n \times n}$.

A key-role in the computational comparison of lattices play the *Hermite and Smith normal forms*. The Hermite normal form defines a canonical choice of primitive vectors, whereas the Smith normal form allows us to find least common sublattices. Using these canonical forms is a crucial ingredient in both the theoretical analysis and the automation of the LFA.

\[def:unimodular\] A matrix $U \in {
    \mathbb{Z}}^{n\times n}$ is called *unimodular* if $\det(U)\in \{\pm 1\}$.

\[def:HNF\] A matrix $H\in {
    \mathbb{R}}^{n\times n}$ is in *Hermite normal form* (HNF) if it is upper triangular, elementwise non-negative and its row-wise maximum is located on the diagonal.

\[def:SNF\] A matrix $S\in {
    \mathbb{R}}^{n\times n}$ is in *Smith normal form* (SNF) if it is a diagonal matrix and the diagonal entries satisfy $ \nicefrac{s_{i+1}}{s_{i}}\in {
    \mathbb{Z}}$ for all $i=1,\ldots,n-1$. These entries are called the *elementary divisors*.

\[thm:integer\_normalforms\] Let ${\mathcal{A}}\in {
    \mathbb{Q}}^{n\times n}$, then

1.  there exists a unimodular $U \in {
        \mathbb{Z}}^{n \times n}$ such that $H = {\mathcal{A}}U$ is in HNF,

2.  there exist unimodular $U,V \in {
        \mathbb{Z}}^{n \times n}$ such that $S = V {\mathcal{A}}U$ is in SNF.

In addition both normal forms $H$ and $S$ are unique with respect to ${\mathcal{A}}$.

Polynomial algorithms to compute the HNF and SNF can for example be found in [@ComputerAlgebraHandbook]. Implementations for the computation of these normalforms are for example part of the PARI software package [@PARI2].

Using these definitions and results one obtains a precise statement about the equality of two lattices.

\[thm:lattice\_equiv\_unimod\] Let ${
    \mathbb{L}}({\mathcal{A}})$ and ${
    \mathbb{L}}({\mathcal{C}})$ be two lattices then the following statements are equivalent.

1.  ${
        \mathbb{L}}({\mathcal{A}}) = {
        \mathbb{L}}({\mathcal{C}})$.

2.  ${\mathcal{A}}^{-1} {\mathcal{C}}\in {
        \mathbb{Z}}^{n\times n}$ and ${\mathcal{C}}^{-1} {\mathcal{A}}\in {
        \mathbb{Z}}^{n\times n}$.

3.  There exists a unimodular matrix $U$, such that ${\mathcal{C}}={\mathcal{A}}U$.

4.  ${\mathcal{A}}$ and ${\mathcal{C}}$ have the same HNF.

Instead of analyzing infinite lattice/crystal structures, we limit ourselves to analyze finite dimensional periodic structures due to the fact that we are ultimately aiming to analyze finite dimensional problems. To this end, another helpful tool is the definition of crystal tori which are defined as quotient groups.

Let ${
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$ be a crystal and ${
    \mathbb{L}}({\mathcal{C}}) \subset {
    \mathbb{L}}({\mathcal{A}})$ be a sublattice. We define the *crystal torus* $T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}$ by $$\begin{aligned}
            T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}:={{\raisebox{.2em}{${
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$}\left/\raisebox{-.2em}{${
    \mathbb{L}}({\mathcal{C}})$}\right.}}.
        \end{aligned}$$ For every $x+{\mathfrak{s}}\in {
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$, their equivalence class $[x+{\mathfrak{s}}]$ is in $T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}$. Furthermore, the elements of $T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}$ are defined by the equivalence $$\begin{aligned}
= [y+{\mathfrak{s}}] \quad \Longleftrightarrow \quad \text{there exists } z \in {
    \mathbb{L}}({\mathcal{C}}), \text{ such that } x = y + z.         \end{aligned}$$

For theoretical and practical reasons, e.g., in \[thm:wavefunction\_ONB\_vec,alg:elements\_in\_quotient\_space\], it is necessary to be able to list all elements of a torus $T^{\mathfrak{s}}_{{\mathcal{A}}, {\mathcal{A}}M } = {\{[x]\in T^{\mathfrak{s}}_{{\mathcal{A}}, {\mathcal{A}}M }  \, : \, x\in {{
    \mathcal{P}}}({\mathcal{A}}M) \cap {
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})\}}$, $M \in {
    \mathbb{{
    \mathbb{Z}}}}^{n \times n}$, uniquely.

To illustrate this point, consider for example an arbitrary lattice ${
    \mathbb{L}}({\mathcal{A}})$ with ${\mathcal{A}}\in \mathbb{R}^{2\times 2}$, and the lattice torus $T_{{\mathcal{A}}, {\mathcal{A}}M }$ with $$M = \mat[c]{m_1 & m_2} =
    \mat[c]{
        2 & 3 \\ 2 & -2
    },\ m_{1},m_{2} \in \mathbb{Z}^{2},$$ as depicted in \[fig:HNF\_LIST\]. Even though we know that the quotient space consists of $|T_{{\mathcal{A}}, {\mathcal{A}}M }|=|\operatorname{det}(M)|=10$ different elements, there is no apparent canonical list of these elements. Fortunately, a canonical ordering of the lattice points on a torus can be formulated using the Hermite or Smith normal form of $M$.

\[cor:quotient\_latticepointlist\] Let $T_{{\mathcal{A}},{\mathcal{C}}}$ be arbitrary, i.e., ${\mathcal{C}}= {\mathcal{A}}M$ for some $M \in {
    \mathbb{{
    \mathbb{Z}}}}^{n \times n}$.

-   Let $H \in{
        \mathbb{{
        \mathbb{Z}}}}^{n \times n}$ be the HNF of $M$ with entries $H_{ij}$ (cf. \[def:HNF\]). Defining the index set $I=I_1 \times I_2 \times \ldots \times I_n$ by $I_\ell:=\{0,1,\ldots,H_{\ell\ell}-1\}$, we then obtain $$\begin{aligned}
                    T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}} & = {\{[x_j+{\mathfrak{s}}] \, : \,  x_j = {\mathcal{A}}j,\ j \in I\}}
                \end{aligned}$$ $\text{with } [x_j+{\mathfrak{s}}] \neq [x_{j'}+{\mathfrak{s}}] \ \Leftrightarrow\ j \neq j' \in I.$

-   Let $S=U^{-1}MV \in {
        \mathbb{Z}}^{n\times n}$ denote the Smith decomposition of $M$ with diagonal entries $s_{i}$ (cf. \[def:SNF\]) and unimodular matrices $U,V$. Defining the index set $\tilde{I}=\tilde{I}_1 \times \tilde{I}_2 \times \ldots \times \tilde{I}_n$ by $\tilde{I}_\ell:=\{0,1,\ldots,s_{\ell}-1\}$, we then obtain $$\begin{aligned}
                    T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}} & = {\{[x_j+{\mathfrak{s}}] \, : \,  x_j = \tilde{{\mathcal{A}}} j,\ j \in \tilde{I}\}}\ 
                \end{aligned}$$ $\text{with }\ [x_j+{\mathfrak{s}}] \neq [x_{j'}+{\mathfrak{s}}] \ \Leftrightarrow\ j \neq j' \in \tilde{I},$ where $\tilde{{\mathcal{A}}} := {\mathcal{A}}U$ denotes the altered lattice basis.

Both statements are a direct consequence of the triangular or diagonal shape of the normal forms and \[thm:lattice\_equiv\_unimod\], i.e., lattices are not changed by unimodular column transformations. On the one hand we have $T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}} = T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{A}}M } = T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{A}}H }.$ The second statement follows from $\tilde{{\mathcal{A}}} S = {\mathcal{C}}V$ and hence $T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}} = T^{\mathfrak{s}}_{{\mathcal{A}}U , {\mathcal{C}}V} = T^{\mathfrak{s}}_{\tilde{{\mathcal{A}}} ,\tilde{{\mathcal{A}}} S }.$

In terms of the example of \[fig:HNF\_LIST\] we obtain:

-   The Hermite normal form $H$ of $M$ is given by $$H= \mat[c]{h_{1} & h_{2}} = 
            \mat[c]{
                5 & 2 \\ 0 & 2},\ h_{1},h_{2} \in \mathbb{Z}^{2}
            .$$ Thus, a unique list of all representatives of $T_{{\mathcal{A}},{\mathcal{A}}M}$ is given by $$\begin{aligned}
                T_{{\mathcal{A}},{\mathcal{A}}M}= T_{{\mathcal{A}},{\mathcal{A}}H} &= {\{[x]=[j_1 {\mathcal{a}}_1 + j_2 {\mathcal{a}}_2] \, : \, (j_1,j_2) \in \{0,1,2,3,4\}\times\{0,1\} \}}. \\
                                  \end{aligned}$$

-   The Smith decomposition of $M$ is given by $$S= \mat[c]{s_{1} & 0 \\ 0 & s_{2}} = 
            \mat[c]{
                1 & 0 \\ 0 & 10}
            =\mat[r]{
                1 & 0 \\ -4 & -1} M \mat[r]{
                -1 & -3 \\ 1 & 2}
            .$$ Thus, another unique list of all representatives of $T_{{\mathcal{A}},{\mathcal{A}}M}$ is given by $$T_{{\mathcal{A}},{\mathcal{A}}M}= T_{\tilde{{\mathcal{A}}} , \tilde{{\mathcal{A}}} S}={\{[x]=[j_1 \tilde{{\mathcal{a}}}_1 + j_2 \tilde{{\mathcal{a}}}_2] \, : \, (j_1,j_2) \in \{0\}\times\{0,1,\ldots,9\} \}},                $$ where $\tilde{{\mathcal{A}}} = \mat{\tilde{{\mathcal{a}}}_1 & \tilde{{\mathcal{a}}}_2} = {\mathcal{A}}\mat[r]{
                1 & 0 \\ -4 & -1}$.

All tori representations $T_{{\mathcal{A}},{\mathcal{A}}M}$, $T_{{\mathcal{A}},{\mathcal{A}}H}$ and $T_{\tilde{{\mathcal{A}}} , \tilde{{\mathcal{A}}} S} $ are depicted in \[fig:HNF\_LIST\]. In the remainder we drop the bracket notation for reasons of readability.

Operators on Crystals {#sec:operators}
=====================

Now that the basic notation of the underlying structure is in place we introduce notation for value distributions and operators on these structures. For aforementioned reasons, we restrict ourselves to the finite dimensional setting by only considering quotient groups of lattices $$\begin{aligned}
        T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}: = {{\raisebox{.2em}{${
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$}\left/\raisebox{-.2em}{${
    \mathbb{L}}({\mathcal{Z}})$}\right.}}
    \end{aligned}$$ with ${
    \mathbb{L}}({\mathcal{Z}}) \subset {
    \mathbb{L}}({\mathcal{A}})$ being an arbitrary sublattice of ${
    \mathbb{L}}({\mathcal{A}})$. While it is significantly easier to be mathematically precise in this general finite setting than in an infinite setting it is also much closer to the targeted applications, namely finite dimensional approximations of PDEs and Hamiltonians. Eventually we only work with operators as part of numerical simulations, i.e., we face only a finite number of unknowns/lattice points anyway. The quotient group $T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}$, ${
    \mathbb{L}}({\mathcal{Z}}) \subset {
    \mathbb{L}}({\mathcal{A}})$, precisely describes such an arbitrarily large but finite torus. To shorten notation we use $T_{{\mathcal{A}}}^{\mathfrak{s}}$ instead of $T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}$ whenever we do not specify ${\mathcal{Z}}$ explicitly.

A crystal-operator is a linear function $$\begin{aligned}
            L : {
    \mathbcal{L}}( T_{{\mathcal{A}}}^{{\mathfrak{d}}} ) \longrightarrow {
    \mathbcal{L}}( T_{{\mathcal{A}}}^{{\mathfrak{c}}} ),
        \end{aligned}$$ where $T_{{\mathcal{A}}}^{{\mathfrak{d}}}$ corresponds to the crystal of the domain and $T_{{\mathcal{A}}}^{{\mathfrak{c}}}$ corresponds to the crystal of the codomain. The function spaces are defined by $$\begin{aligned}
            {
    \mathbcal{L}}( T_{{\mathcal{A}}}^{{{{\mathfrak{s}}}}} ) = {\{f=(f_1,\ldots,f_{|{{\mathfrak{s}}}|})  \, : \, T_{{\mathcal{A}}}\longrightarrow {
    \mathbb{C}}^{|{{{{\mathfrak{s}}}}}|}\}}.
        \end{aligned}$$ A value $f_j(x)$, $x\in{
    \mathbb{L}}({\mathcal{A}})$, corresponds to a value at the position $x + {\mathfrak{s}}_j$. The function space is equipped with the scalar product $${\langle f,g\rangle_{}} := \frac{1}{|T_{{\mathcal{A}},{\mathcal{Z}}}|} \sum_{x\in T_{{\mathcal{A}},{\mathcal{Z}}}} {\langle f(x),g(x)\rangle_{2}},$$ where ${\langle f(x),g(x)\rangle_{2}} := \sum_{\ell=1}^{|{{\mathfrak{s}}}|} f_\ell(x) \overline{g_\ell(x)}$ denotes the Euclidean scalar product on ${
    \mathbb{C}}^{|{{\mathfrak{s}}}|}$.

In the context of LFA we are interested in operators which can be represented in (block) stencil notation. That is, translationally invariant operators that can be written as multiplication operators. As these two properties are in fact equivalent (cf. [@nla.cat-vn1897929 Theorem 3.16]), we can connect those operators to the notation of crystal structures.

Let $L : {
    \mathbcal{L}}( T_{{\mathcal{A}}}^{{\mathfrak{d}}} ) \longrightarrow {
    \mathbcal{L}}( T_{{\mathcal{A}}}^{{\mathfrak{c}}} )$ be a crystal operator. The following statements are equivalent.

1.  $L$ is a multiplication operator, i.e., there exist matrices $m_L^{(y)} \in {
        \mathbb{C}}^{|{{\mathfrak{c}}}|\times |{{\mathfrak{d}}}|}$ such that for each $x\in T_{\mathcal{A}}$ and $f\in {
        \mathbcal{L}}( T_{{\mathcal{A}}}^{{\mathfrak{d}}} ) $ we have $$\begin{aligned}
                    (Lf)(x) = \sum_{y \in T_{\mathcal{A}}} m_L^{(y)} f(x+y).
                \end{aligned}$$

2.  $L$ is $({\mathcal{A}})$-translationally invariant, i.e., $${L}{\mathbb{T}}_{\mathcal{a}}- {\mathbb{T}}_{\mathcal{a}}{L}= 0 \quad \text{for all (primitive) vectors } {\mathcal{a}}\in {
        \mathbb{L}}({\mathcal{A}}),$$ where the translation operator is defined by $({\mathbb{T}}_{\mathcal{a}}f)(x) = f(x+{\mathcal{a}})$.

For the analysis of such operators the concept of the dual lattice comes in handy as already considered in similar form in [@GaspGracLisb2009; @ZhouFult].

\[def:dual\_lattice\] Let ${
    \mathbb{L}}({\mathcal{A}})$ be a lattice. Its *dual lattice* ${
    \mathbb{L}}({\mathcal{B}})={
    \mathbb{L}}({\mathcal{A}})^*$ is the set $$\begin{aligned}
            {
    \mathbb{L}}({\mathcal{A}})^*:={\{k \in {
    \mathbb{R}}^n \, : \,  {\langle k,x\rangle_{2}} \in {
    \mathbb{Z}}\text{ for all } x\in{
    \mathbb{L}}\}}.
        \end{aligned}$$ A lattice basis of the dual lattice is given by ${\mathcal{B}}= {\mathcal{A}}^{-T}$. The elements of ${
    \mathbb{L}}({\mathcal{A}})^*$ may also be referred to as *wave vectors*.

In addition to the dual space we introduce an orthonormal basis of wave functions that are compatible with the crystal structure introduced in \[sec:lattices\_crystals\]. This basis is an extension of the basis used in [@Kuo1989] to analyze the red-black Gauss-Seidel relaxation. Furthermore, similar basis functions have been used in the context of LFA, e.g. in [@BoonLentVand2008; @1803.08864; @RittichBolten2018; @1902.10248; @Brown2018LOCALFA].

\[thm:wavefunction\_ONB\_vec\] An orthonormal basis for the function space ${
    \mathbcal{L}}(T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{Z}}})$ with a structure element ${\mathfrak{s}}=({\mathfrak{s}}_1,\ldots,{\mathfrak{s}}_m)$ is given by the wave functions $e_{1,k},e_{2,k},\ldots,e_{m,k}$ defined by $$\begin{aligned}
            \left(e_{\ell,k}(x)\right)_j := 
            \begin{cases}
                {e^{2\pi i{\langle k,x\rangle_{2}}}} & \text{ if }j=\ell,                 \\ 0 & \text{ else,} 
            \end{cases}  
        \end{aligned}$$ with $k \in T_{{\mathcal{A}},{\mathcal{Z}}}^*:= {{\raisebox{.2em}{${
    \mathbb{L}}({\mathcal{Z}})^*$}\left/\raisebox{-.2em}{${
    \mathbb{L}}({\mathcal{A}})^*$}\right.}} $.

The statement follows by a straightforward, but lengthy calculation, by assuming without loss of generality that ${\mathcal{A}}^{-1}{\mathcal{Z}}$ is given in Smith normal form, making use of \[cor:quotient\_latticepointlist\] and the geometric sum formula.

An illustration of a lattice torus $T_{{\mathcal{A}},{\mathcal{Z}}}$ along with its dual $T_{{\mathcal{A}},{\mathcal{Z}}}^*$ is given in \[fig:TorusAndDual\].

The orthonormal basis of \[thm:wavefunction\_ONB\_vec\] can be split into subsets with respect to the wave vector $k$, i.e., ${
    \mathbcal{L}}(T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{Z}}}) = \cup_{k \in T_{{\mathcal{A}},{\mathcal{Z}}}^*} \operatorname{span}(H_k) $ with $$\label{eq:space-of-harmonics}
        H_k = \operatorname{span}{\{e_{\ell,k} \, : \, \ell=1,\ldots,m\}}.$$

\[thm:L-invariance\] Let ${L}: {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}})$ be a multiplication operator. Then the subspaces $H_k$ of \[eq:space-of-harmonics\] are $L$-invariant, i.e., $L(H_k) \subseteq H_k$.

Let $e_{k}(x)$ denote an arbitrary value distribution in $H_{k}$. That is, there exist $\alpha_{1},\ldots,\alpha_{m} \in {
    \mathbb{C}}$ such that $$\begin{aligned}
            e_k(x)=\sum_{\ell}\alpha_\ell e_{\ell,k} = 
                        (\alpha_1 {e^{2\pi i{\langle k,x\rangle_{2}}}},\ldots,\alpha_m{e^{2\pi i{\langle k,x\rangle_{2}}}})^T \in \operatorname{span}(H_k).
        \end{aligned}$$ Then we obtain by direct calculation $$\label{eq:symbol}
            (L e_k)(x)  = \sum_{y} m_L^{(y)} e_k(x+y) = \left(\sum_{y} m_L^{(y)} {e^{2\pi i{\langle k,y\rangle_{2}}}}\right) e_k(x).$$

Due to their $L$-invariance the subspaces $H_k$ are oftentimes referred to as *spaces of harmonics*. Thus we can easily represent any ${\mathcal{A}}$-translationally operator via its *symbols*, which are formally defined as follows.

\[def:crystal\_op\_symbol\] Let $L:{
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{t}})$ be a multiplication operator with $$\begin{aligned}
                                    (Lf)(x) = \displaystyle\sum_{y \in T_{{\mathcal{A}},{\mathcal{Z}}}} m_L^{(y)} f(x+y),\  m_L^{(y)} \in {
    \mathbb{C}}^{|{\mathfrak{t}}|\times |{\mathfrak{s}}|}.
                    \end{aligned}$$ We define the *symbol* of $L$ according to \[eq:symbol\] by $$\begin{aligned}
            L_k:=\sum_{y \in T_{{\mathcal{A}},{\mathcal{Z}}}} m_L^{(y)} m_k^{(y)}
                        \text{ with }
                        m_k^{(y)} :=  {e^{2\pi i{\langle k,y\rangle_{2}}}}.                     \end{aligned}$$

In case ${\mathfrak{s}}= {\mathfrak{t}}$ the spectrum of $L$ can then be extracted from its symbols $L_{k}$.

\[cor:eigenvalues\_crystal\_operator\] Let $L:{
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}})$ be a multiplication operator with $$\begin{aligned}
                                                (Lf)(x) = \sum_{y \in T_{{\mathcal{A}},{\mathcal{Z}}}} m_L^{(y)} f(x+y),\ m_L^{(y)} \in {
    \mathbb{C}}^{|{\mathfrak{s}}|\times |{\mathfrak{s}}|}.
                    \end{aligned}$$ Then ${\operatorname{spec}}(L) = \cup_{k\in T_{{\mathcal{A}},{\mathcal{Z}}}^*} {\operatorname{spec}}(L_k)$.

Follows immediately due to the orthonormality of the basis $e_{\ell,k}$ (cf. \[thm:wavefunction\_ONB\_vec\]) and the $L$-invariance of the subspaces $H_{k}$ (cf. \[thm:L-invariance\]).

\[rem:freq\_sampling\] The main purpose of ${\mathcal{Z}}$, i.e., the set of primitive vectors that define an arbitrary sublattice ${
    \mathbb{L}}({\mathcal{Z}})$ of ${
    \mathbb{L}}({\mathcal{A}})$, is to simplify the theory developed in \[sec:operators\] by turning an infinite dimensional setting to an (arbitrarily large) finite one. Though there is another interpretation to it as well. In \[cor:eigenvalues\_crystal\_operator\] ${\mathcal{Z}}$ explicitly specifies the resolution of the frequency space as seen in \[fig:TorusAndDual\], i.e., $$T_{{\mathcal{A}},{\mathcal{Z}}}^* = {
    \mathbb{L}}({\mathcal{Z}}^{-T}) \cap {{
    \mathcal{P}}}({\mathcal{A}}^{-T}),$$ where the spectrum of the multiplication operator is sampled. Due to the reciprocal nature of the dual space, the larger $|\operatorname{det}({\mathcal{Z}})|$ is, the finer the resolution becomes.

With these tools at hand we are able to fully analyze a single multiplication operator, but typically we are interested in analyzing a composition of several operators using LFA. As long as the corresponding domains and codomains of these operators are compatible we can use the rules of computation given in \[sec:rules\_of\_comp\] on the level of multiplication operators and/or on the level of the corresponding symbols. While computing the sum, a product and taking the transpose can easily be done on both levels, taking the (pseudo-)inverse is simple only on the level of symbols. The (pseudo-)inverse of a multiplication operator may have an arbitrarily large number[^6] of multipliers $m_{L^{-1}}^{(y)} \neq 0$ and thus there is no simple rule to compute it. In case a pseudo-inverse has to be used in the following we opt to employ the Moore-Penrose pseudo-inverse, which we denote by $S^{\dagger}$ for a multiplication operator $S$.

In our framework, which tries to automate as much of the LFA as possible, we would like to allow for user friendly descriptions of all occurring operators. That is, it should be possible to describe operators in terms of their individual translational invariance and ordering of the structure element without having to worry about compatibility issues with other operators on the input level of the analysis. The process how to automatically make crystal representations of operators compatible is explained in detail in \[sec:crystal\_reps\_and\_isos\]. Before diving into the gritty details of this automation process we would like to illustrate the developments made so far with an example in order to convey the introduced notation.

An example {#sec:example}
==========

Consider the red-black Gauss-Seidel method applied to the discretized Laplacian on the unit square with periodic boundary conditions. The most fundamental representation of the discretized unit square with periodic boundary conditions is given by the torus $T_{{\mathcal{A}},{\mathcal{Z}}}$ with $$\begin{aligned}
        {\mathcal{Z}}=\mat[c]{1 & 0 \\ 0 & 1}, {\mathcal{A}}=\mat[c]{{\mathcal{a}}_1 & {\mathcal{a}}_2}=\frac{1}{h} \mat[c]{1 & 0 \\ 0 & 1}\text{ for some }h=\frac{1}{N}, N\in{
    \mathbb{N}}.
    \end{aligned}$$ Then, the discretized Laplacian $L_h : {
    \mathbcal{L}}( T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)} ) \longrightarrow {
    \mathbcal{L}}( T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)} )$ using finite differences is given by $(L_h f)(x) := \sum_{y\in{
    \mathbb{L}}({\mathcal{A}})} m_{L_h}^{(y)} f(x+y)$ with non-zero multipliers:

\(c) at (0,0) [$=$]{}; ; ;

\(c) at (0,) [$=$]{}; ; ;

\(c) at (0,-) [$=$]{}; ; ;

\(c) at (2.2\*,0) [$=$]{}; ; ;

\(c) at (-1.9\*,0) [$=$]{}; ; ;

The error propagator $G$ of the red-black Gauss-Seidel method can be written with respect to the crystal $T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)}$ via $$\begin{aligned}
        G=(I-S_b^\dagger L_h)(I-S_r^\dagger L_h) : {
    \mathbcal{L}}( T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)} ) \longrightarrow {
    \mathbcal{L}}( T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)} )
    \end{aligned}$$ with $$\begin{aligned}
        (S_r f)(x) :=
        \begin{cases}
            \frac{4}{h^2} f(x) & x\in X_\text{red}   \\
            0                  & x\in X_\text{black}
        \end{cases}
        \text{ and } (S_b f)(x) :=
        \begin{cases}
            0                  & x\in X_\text{red}   \\
            \frac{4}{h^2} f(x) & x\in X_\text{black}
        \end{cases}
        ,
    \end{aligned}$$ where $X_\text{red}$ and $X_\text{black}$ correspond to the red and black unknowns of the torus $T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)}$ as illustrated in \[fig:redblack\_laplacian\], b). In order to analyze this composite operator in our framework, we write all occurring operators as multiplication operators. It is now important that the red-black splitting of $T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)} = X_{\text{red}} \cup X_{\text{black}}$ implies the crystal $ T_{{\mathcal{C}},{\mathcal{Z}}}^{\mathfrak{s}}$ with $$\begin{aligned}
        {\mathcal{c}}_1 = {\mathcal{a}}_1 + {\mathcal{a}}_2, \text{ } {\mathcal{c}}_2 = {\mathcal{a}}_1 - {\mathcal{a}}_2, \ {\mathfrak{s}}:= (0, {\mathcal{a}}_1)  = {
    \mathbb{L}}({\mathcal{A}}) \cap {{
    \mathcal{P}}}({\mathcal{C}})
    \end{aligned}$$ such that $T_{{\mathcal{C}},{\mathcal{Z}}}^{(0)} =  X_{\text{red}}$, $T_{{\mathcal{C}},{\mathcal{Z}}}^{({\mathcal{a}}_1)} = X_{\text{black}}$ and $T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)} = X_{\text{red}} \cup X_{\text{black}} \cong T_{{\mathcal{C}},{\mathcal{Z}}}^{\mathfrak{s}}$ (cf. \[fig:redblack\_laplacian\], b)). With respect to this crystal the operators $S_r$ and $S_b$ can be written as multiplication operators $\hat{S}_r, \hat{S}_b : {
    \mathbcal{L}}( T_{{\mathcal{C}},{\mathcal{Z}}}^{{\mathfrak{s}}} ) \longrightarrow {
    \mathbcal{L}}( T_{{\mathcal{C}},{\mathcal{Z}}}^{{\mathfrak{s}}} )$ with $$\begin{aligned}
        (\hat{S}_r f)(x) = \mat[c]{\frac{4}{h^2} & 0 \\ 0 & 0} f(x) \text{ and } (\hat{S}_b f)(x) = \mat[c]{0 & 0 \\ 0 & \frac{4}{h^2}} f(x).
    \end{aligned}$$

We now have described all individual operators of $G$, each one defined with respect to its own (minimal) translational invariance, but the domains and codomains are not identical, i.e., ${\mathcal{A}}\neq {\mathcal{C}}$, such that we cannot directly use the computation rules described in \[sec:rules\_of\_comp\]. In order to carry on with the analysis we have to rewrite the operators with respect to a common crystal structure. In this example we construct this structure by hand, but this process can be automated as explained in \[sec:crystal\_reps\_and\_isos\].

As the crystal $T_{{\mathcal{C}},{\mathcal{Z}}}^{\mathfrak{s}}$ is yet another representation of $T_{{\mathcal{A}},{\mathcal{Z}}}^{(0)}$, the operator $L_h$ can be rewritten with respect to this crystal (cf. \[fig:redblack\_laplacian\]) as $ \hat{L}_h: {
    \mathbcal{L}}( T_{{\mathcal{C}},{\mathcal{Z}}}^{{\mathfrak{s}}} ) \longrightarrow {
    \mathbcal{L}}( T_{{\mathcal{C}},{\mathcal{Z}}}^{{\mathfrak{s}}} )$ with $(\hat{L}_h f)(x) := \sum_{y\in{
    \mathbb{L}}({\mathcal{C}})} m_{\hat{L}_h}^{(y)} f(x+y)$ and non-zero multipliers:

Thus, the spectrum of the error propagator of the red-black Gauss-Seidel method applied to the Laplacian $$\label{eq:RBGEP}
        \hat{G} = (I-\hat{S}_b^\dagger \hat{L}_h)(I-\hat{S}_r^\dagger \hat{L}_h)$$ can now be obtained elementwise for each fixed $k \in T_{{\mathcal{C}},{\mathcal{Z}}}^* = {{
    \mathcal{P}}}({\mathcal{C}}^{-T}) \cap {
    \mathbb{L}}({\mathcal{Z}}^{-T})$ by first computing the individual *symbols* $I_k, (\hat{S}_r)_k, (\hat{S}_b)_k$ and $(\hat{L}_h)_k$ and assembling the symbols $\hat{G}_k$ according to \[eq:RBGEP\] and the rules in \[sec:rules\_of\_comp\], followed by the computation of the eigenvalues of the matrices $\hat{G}_k$. The resulting spectral information of the discretized Laplace operator $\hat{L}_h$ and the error propagator $\hat{G}$ is illustrated in \[fig:redblack\_spectrum\], where it is sampled on the dual lattice $T_{{\mathcal{C}},{\mathcal{Z}}}^*$. Note, that one naturally obtains two eigenvalues per sampled wave vector $k$. In case of the spectrum of the red-black Gauss-Seidel error propagator, one of the two eigenvalues is equal to zero for all wave vectors $k$.

Crystal representations and natural isomorphisms {#sec:crystal_reps_and_isos}
================================================

In general, we are given several multiplication operators which make up the error propagator of an iterative method, each defined with respect to its own (minimal) translational invariance. In order to analyze the method we thus need to find a common denominator, i.e., a lattice basis corresponding to the collective translational invariance, and rewrite the operators accordingly. The following theorem yields a set of primitive vectors of such a collective translational invariance for two arbitrary lattices, if it exists.

\[thm:lcm\_mat\] Given two $n$-dimensional lattices ${
    \mathbb{L}}({\mathcal{A}})$, ${
    \mathbb{L}}({\mathcal{B}})$. If there exists an $r\in {
    \mathbb{Z}}$, such that $M = r {\mathcal{A}}^{-1} {\mathcal{B}}\in {
    \mathbb{Z}}^{n\times n}$, then there is a lattice ${
    \mathbb{L}}({\mathcal{C}})$ with ${
    \mathbb{L}}({\mathcal{C}}) \subset {
    \mathbb{L}}({\mathcal{A}})$ and ${
    \mathbb{L}}({\mathcal{C}}) \subset {
    \mathbb{L}}({\mathcal{B}})$ with $|\det({\mathcal{C}})|$ as small as possible. A lattice basis of ${
    \mathbb{L}}({\mathcal{C}})$ is given by $$\begin{aligned}
            {\mathcal{C}}= {\mathcal{B}}T^{-1} N_{\mathcal{B}}, \quad
        \end{aligned}$$ where $N_{\mathcal{B}}$ is a diagonal matrix with $(N_{\mathcal{B}})_{i,i} := r \cdot \operatorname{gcd}(r,s_{i})^{-1}$, where $S= V^{-1} M T^{-1} = \operatorname{diag}{(s_{1},\ldots,s_{n})}$ is the SNF of $M$ (cf. \[def:SNF\]) and $\operatorname{gcd}(r,s_{i})$ denotes the greatest common divisor of $r$ and $s_i$. Consequently, we write ${
    \mathbb{L}}({\mathcal{C}}) = \operatorname{lcm}( {
    \mathbb{L}}({\mathcal{A}}), {
    \mathbb{L}}({\mathcal{B}}))$ and call it the *least common multiple* of ${
    \mathbb{L}}({\mathcal{A}})$ and ${
    \mathbb{L}}({\mathcal{B}})$.

Due to \[lem:sublattice\] it is sufficient to find integral matrices $N_{\mathcal{A}}, N_{\mathcal{B}}$, such that $$\begin{aligned}
                {
    \mathbb{L}}({\mathcal{C}}) = {
    \mathbb{L}}({\mathcal{A}}N_{\mathcal{A}}) = {
    \mathbb{L}}({\mathcal{B}}N_{\mathcal{B}})
            \end{aligned}$$ with $|\det(N_{\mathcal{A}})|$ and $|\det(N_{\mathcal{B}})|$ as small as possible. Using \[thm:lattice\_equiv\_unimod\], i.e., $
            {
    \mathbb{L}}({\mathcal{B}}) = {
    \mathbb{L}}({\mathcal{B}}V_1), \quad {
    \mathbb{L}}({\mathcal{C}}) = {
    \mathbb{L}}({\mathcal{C}}U_1),
            $ for any unimodular matrices $U_1,V_1$, we can assume the equality $$\begin{aligned}
                {\mathcal{A}}N_{\mathcal{A}}= {\mathcal{B}}U N_{\mathcal{B}}\end{aligned}$$ for any unimodular $U$ and $N_{\mathcal{B}}$ in Hermite normal form (cf. \[thm:integer\_normalforms\]). Plugging in the Smith decomposition $VST$ of $M=r {\mathcal{A}}^{-1} {\mathcal{B}}$ and defining $U:=T^{-1}$, we find $$\begin{aligned}
                \begin{array}{rcccl}
                    N_{\mathcal{A}}& = & {\mathcal{A}}^{-1} {\mathcal{B}}T^{-1} N_{\mathcal{B}}&= &  \frac{1}{r} VS N_{\mathcal{B}}.     
                \end{array}
            \end{aligned}$$ Both matrices $$\begin{aligned}
                N_{\mathcal{B}}= \mat{(N_{\mathcal{B}})_{1,1} & \ldots & (N_{\mathcal{B}})_{1,n} \\
                    & \ddots & \vdots \\
                    & & (N_{\mathcal{B}})_{n,n} } \text{ and } \frac{1}{r} S N_{\mathcal{B}}= \mat{\frac{s_{1}}{r}(N_{\mathcal{B}})_{1,1} & \ldots & (N_{\mathcal{B}})_{1,n}\frac{s_{1}}{r} \\
                    & \ddots & \vdots \\
                    & & (N_{\mathcal{B}})_{n,n}\frac{s_{n}}{r} }
            \end{aligned}$$ have to be integral with $$|\det(N_{\mathcal{B}})| = |\prod_{i=1}^n (N_{\mathcal{B}})_{i,i}| \text{ and }|\det(\frac{1}{r} S N_{\mathcal{B}})| = |\prod_{i=1}^n (N_{\mathcal{B}})_{i,i} \frac{s_{i}}{r} |$$ as small as possible. It can easily be verified that $$\begin{aligned}
                (N_{\mathcal{B}})_{i,i} := \frac{r}{\operatorname{gcd}(r,s_{i})}
            \end{aligned}$$ is the optimal choice for the diagonal entries. With this choice, the off-diagonal entries $(N_{\mathcal{B}})_{i,j}$ have to be integral multiples of $(N_{\mathcal{B}})_{i,i}$. Due to the fact that $N_{\mathcal{B}}$ is in Hermite normal form, the off-diagonal entries are zero.

We now study different representations of the same crystal structure in order to derive a general way to rewrite a multiplication operator with respect to some coarser crystal structure corresponding to a sublattice, as has been done manually in \[sec:example\] for the discretized Laplacian.

\[thm:crystal\_congruence\] Let ${
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$ be a crystal and ${
    \mathbb{L}}({\mathcal{C}}) \subset {
    \mathbb{L}}({\mathcal{A}})$ a sublattice. Denoting $T_{{\mathcal{A}},{\mathcal{C}}}=\{{\mathfrak{t}}_1,\ldots,{\mathfrak{t}}_p\}$, [^7] the set $$\begin{aligned}
            {\Xi}({\mathcal{C}}):={\{x+\delta\in{
    \mathbb{R}}^n \, : \, \delta\in{{
    \mathcal{P}}}({\mathcal{A}}), x \in \{{\mathfrak{t}}_1,\ldots,{\mathfrak{t}}_p\}\}}
        \end{aligned}$$ defines a primitive cell of ${
    \mathbb{L}}({\mathcal{C}})$, and the tuple $$\begin{aligned}
            {\mathfrak{u}}= ({\mathfrak{t}}_1+{\mathfrak{s}}_1,\ldots , {\mathfrak{t}}_1+{\mathfrak{s}}_m, {\mathfrak{t}}_2+{\mathfrak{s}}_1,\ldots,{\mathfrak{t}}_p+{\mathfrak{s}}_m)\in{\Xi}({\mathcal{C}})^{p\cdot m}
        \end{aligned}$$ defines a structure element of ${
    \mathbb{L}}^{\mathfrak{u}}({\mathcal{C}})$ such that ${
    \mathbb{L}}^{\mathfrak{u}}({\mathcal{C}}) \cong {
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}}),$ meant as a one-to-$p$ correspondence.

Without loss of generality we may assume ${\mathfrak{t}}_{j} \in {{
    \mathcal{P}}}({\mathcal{C}})$, $j = 1,\ldots,p$. Then, each element in ${
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$ can be written as $$\begin{aligned}
            z=({\mathcal{A}}x+{\mathfrak{s}}_1,\ldots,{\mathcal{A}}x+{\mathfrak{s}}_m)
        \end{aligned}$$ and there is a unique $y$, such that ${\mathcal{A}}x={\mathcal{C}}y={\mathcal{C}}\lfloor y \rfloor + {\mathcal{C}}(y - \lfloor y \rfloor)$ with ${\mathcal{C}}\lfloor y \rfloor \in {
    \mathbb{L}}({\mathcal{C}})$ and ${\mathcal{C}}(y - \lfloor y \rfloor)={\mathfrak{t}}_j \in {{
    \mathcal{P}}}({\mathcal{C}})\cap{
    \mathbb{L}}({\mathcal{A}})$. Thus, we find $z$ as a unique part of the element $$\begin{aligned}
            {\mathcal{C}}\lfloor y \rfloor +{{\mathfrak{u}}}=(\ldots,{\mathcal{C}}\lfloor y \rfloor  + {\mathfrak{t}}_j + {\mathfrak{s}},\ldots)=(\ldots,z,\ldots).
        \end{aligned}$$ This argument works in the other direction in the same way.

\[rem:natural\_isomorph\] Note, that ${\mathfrak{u}}$, as defined in \[thm:crystal\_congruence\], is an explicit representation of $T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}$, thus using this explicit representation \[thm:crystal\_congruence\] implies a congruence of the function spaces $$\begin{aligned}
        {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}) \cong {
    \mathbcal{L}}(T_{{\mathcal{C}},{\mathcal{Z}}}^{T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}})
    \end{aligned}$$ given by the natural isomorphism $\eta: {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{C}},{\mathcal{Z}}}^{T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}})$, $$\begin{aligned}
        {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{\mathfrak{s}}) \ni f(\cdot) \mapsto (\eta f )(\cdot)= (f(\cdot+{\mathfrak{t}}_1),\ldots, f(\cdot+{\mathfrak{t}}_p)) \in {
    \mathbcal{L}}(T_{{\mathcal{C}},{\mathcal{Z}}}^{T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{C}}}}),
    \end{aligned}$$ as $(f)$ and $(\eta f)$ describe the same value distribution on the crystal. This congruence in turn implies that the coarsest possible crystal interpretation, i.e., ${\mathfrak{u}}\cong T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{Z}}}$, is simply the complex coordinate space $${
    \mathbb{C}}^{mn} = {
    \mathbcal{L}}(T_{{\mathcal{Z}},{\mathcal{Z}}}^{T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{Z}}}})  \cong {
    \mathbcal{L}}(T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{Z}}})$$ with $n:=|T_{{\mathcal{A}},{\mathcal{Z}}}|$. The scalar product on ${
    \mathbcal{L}}(T^{\mathfrak{s}}_{{\mathcal{A}},{\mathcal{Z}}})$ then corresponds to the Euclidean scalar product on ${
    \mathbb{C}}^{mn}$ up to a factor of $n$.

Using the natural isomorphism of function spaces in \[rem:natural\_isomorph\], we can derive the transformations of multiplication operators when coarsening the underlying crystal representation corresponding to a sublattice.

\[thm:operator\_sublattice\_coarsening\] Consider crystals ${
    \mathbb{L}}^{\mathfrak{d}}({\mathcal{A}})$, ${
    \mathbb{L}}^{\mathfrak{c}}({\mathcal{A}})$, a sublattice ${
    \mathbb{L}}({\mathcal{C}}) \subset {
    \mathbb{L}}({\mathcal{A}})$ and a multiplication operator $$\begin{aligned}
            L: {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{d}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{c}}}), \quad ({L}f)(x):=\sum_{y\in T_{{\mathcal{A}}}} m_L^{(y)} f(x+y), \quad m_L^{(y)} \in {
    \mathbb{C}}^{|{\mathfrak{c}}|\times |{\mathfrak{d}}|}.
        \end{aligned}$$ Then, using $T_{{\mathcal{A}},{\mathcal{C}}}=\{{\mathfrak{t}}_1,\ldots,{\mathfrak{t}}_p\}$, the multiplication operator $$\begin{aligned}
            G:{
    \mathbcal{L}}(T_{\mathcal{C}}^{\hat{{\mathfrak{d}}}}) \rightarrow  {
    \mathbcal{L}}(T_{\mathcal{C}}^{\hat{{\mathfrak{c}}}}), \quad (G g)(x)=\sum_{y\in T_{\mathcal{C}}} m_G^{(y)} g(x+y), \quad m_G^{(y)} \in {
    \mathbb{C}}^{p |{\mathfrak{c}}|\times p|{\mathfrak{d}}|},
        \end{aligned}$$ with block matrices $(m_G^{(y)})_{i,k} := m_L^{(y-{\mathfrak{t}}_i+{\mathfrak{t}}_k)} \in {
    \mathbb{C}}^{|{\mathfrak{c}}|\times |{\mathfrak{d}}|}$ fulfills the commutative diagram: $$\begin{tikzcd}
            {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{d}}}) \arrow[r, "{L}"] \arrow[d, "\eta^{\mathfrak{d}}"]
            & {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{c}}}) \arrow[d, "\eta^{\mathfrak{c}}" ] \\
            {
    \mathbcal{L}}(T_{\mathcal{C}}^{\hat{{\mathfrak{d}}}}) \arrow[r, "G" ]
            & {
    \mathbcal{L}}(T_{\mathcal{C}}^{\hat{{\mathfrak{c}}}}).
        \end{tikzcd}$$ Here, the mappings for ${\mathfrak{s}}\in \{{\mathfrak{d}},{\mathfrak{c}}\}$, $$\begin{aligned}
            \eta^{\mathfrak{s}}: {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{s}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{C}}^{\hat{{\mathfrak{s}}}}), \quad f(\cdot)\mapsto (f(\cdot +{\mathfrak{t}}_1),\ldots,f(\cdot+{\mathfrak{t}}_p)),
        \end{aligned}$$ denote the natural isomorphisms between the congruent crystal representations.

A straightforward calculation for each block-row $i$ yields $$\begin{aligned}
            \begin{array}{rcl}
                [(\eta^{\mathfrak{c}}L f)(x)]_i & = & (Lf)(x+{\mathfrak{t}}_i)                                                                                            \ =\ \displaystyle\sum\limits_{k=1}^{p} \displaystyle\sum\limits_{y \in T_{\mathcal{C}}} m_L^{(y+{\mathfrak{t}}_k)} f(x+y+{\mathfrak{t}}_i+{\mathfrak{t}}_k) \\[\medskipamount]                & = &  \displaystyle\sum\limits_{k=1}^{p} \displaystyle\sum\limits_{y \in T_{\mathcal{C}}} m_L^{(y-{\mathfrak{t}}_i+{\mathfrak{t}}_k)} f(x+y+{\mathfrak{t}}_k)                 \ =\   \displaystyle\sum\limits_{y \in T_{\mathcal{C}}} \displaystyle\sum\limits_{k=1}^{p} (m_G^{(y)})_{i,k} f(x+y+{\mathfrak{t}}_k)      \\[\bigskipamount]                & = & [G (f(x+{\mathfrak{t}}_1),\ldots,f(x+{\mathfrak{t}}_p))]_i                                                    \ =\ [(G \eta^{\mathfrak{d}}f)(x)]_i.
            \end{array}
        \end{aligned}$$

Using \[thm:operator\_sublattice\_coarsening\] we now know how to rewrite multiple multiplication operators with respect to some common crystal structure with a coarser translational invariance. Due to the fact that we do not make any assumption on the initial representation of the crystal structures, the resulting structure elements of \[cor:quotient\_latticepointlist\] might differ in their orderings and might contain shifts with respect to the common shift invariance. To automatically remove these differences and determine the corresponding transformations of the associated multiplication operators we first define the notion of congruent structure elements.

\[def:crystaltori\_congruent\] Two crystal tori $T_{\mathcal{A}}^{\mathfrak{s}}\cong T_{\mathcal{A}}^{\mathfrak{t}}$, ${\mathcal{A}}\in {
    \mathbb{R}}^n$ are *congruent* with respect to ${
    \mathbb{L}}({\mathcal{A}})$ if the structure elements are of the same size, i.e., $|{\mathfrak{s}}| = |{\mathfrak{t}}| = m$, and there is a permutation $\pi:\{1,\ldots,m\} \rightarrow \{1,\ldots,m\}$ as well as shifts $y_j \in {
    \mathbb{L}}({\mathcal{A}})$, such that $$\begin{aligned}
            {\mathfrak{s}}_j =  y_j + {\mathfrak{t}}_{\pi(j)}         \end{aligned}$$

In order to introduce a unique representation for the sake of automation, we introduce the following *normal form* and the required transformations to transfer any operator to this form.

\[def:crystal\_op\_and\_normalform\] Let $L: {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{d}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{c}}})$ be a multiplication operator. We say $L$ is in *normal form* if

-   the coordinates of the structure elements are found in the primitive cell spanned by the primitive vectors, i.e., ${\mathfrak{d}}_i,{\mathfrak{c}}_j \in {{
        \mathcal{P}}}({\mathcal{A}}) = {\mathcal{A}}[0,1)^n \text{\ for each\ } i,j,$

-   the structure elements ${\mathfrak{d}}$ and ${\mathfrak{c}}$ are sorted lexicographically.[^8]

We now derive the implications of \[def:crystaltori\_congruent\] for multiplication operators when the structure element is element wise shifted or permuted. We do so in two steps, \[thm:shifted\_structure\_element\] and \[thm:permuted\_structure\_element\]. First, we show that a shift of an entry of the structure element in the codomain or domain results in a modification of the corresponding row or column of the non-zero multipliers, respectively.

\[thm:shifted\_structure\_element\] Consider the two multiplication operators $ L: {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{s}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{t}}})$ and $G: {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{t}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{u}}})$ defined by $$\begin{aligned}
            (Lf)(x) = \sum_{y \in T_{{\mathcal{A}}}} m_L^{(y)} f(x+y), \quad m_L^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{t}}|\times |{\mathfrak{s}}|},  \\
            (Gg)(x) = \sum_{y \in T_{{\mathcal{A}}}} m_G^{(y)} g(x+y), \quad m_G^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{u}}| \times |{\mathfrak{t}}|}.
        \end{aligned}$$ Further let $\hat{{\mathfrak{t}}}$ be a structure element which is obtained from ${\mathfrak{t}}$ when shifted element-wise along ${
    \mathbb{L}}({\mathcal{A}})$, i.e., $$\begin{aligned}
                        {\mathfrak{t}}& =  ({\mathfrak{t}}_1,\ldots,{\mathfrak{t}}_m)    =  (\hat{{\mathfrak{t}}}_1 + y_1,\ldots,\hat{{\mathfrak{t}}}_m+y_m ) = \hat{{\mathfrak{t}}} +(y_1,\ldots,y_m),
        \end{aligned}$$ where $y_1,\ldots,y_m \in{
    \mathbb{L}}({\mathcal{A}})$ and $m=|{\mathfrak{t}}|$. Then, the operators $\hat{L}$ and $\hat{G}$ given by $$\begin{aligned}
            \begin{array}{rcllcl}
                (\hat{L}f)(x) &=& \displaystyle\sum\limits_{y \in T_{{\mathcal{A}}}} m_{\hat{L}}^{(y)} f(x+y), & m_{\hat{L}}^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{t}}|\times |{\mathfrak{s}}|} &\text{ with }& (m_{\hat{L}}^{(y)})_{i,j} := (m_{L}^{(y + y_i)})_{i,j},\\
                (\hat{G}f)(x) &=& \displaystyle\sum\limits_{y \in T_{{\mathcal{A}}}} m_{\hat{G}}^{(y)} f(x+y), & m_{\hat{G}}^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{u}}|\times |{\mathfrak{t}}|} &\text{ with }& (m_{\hat{G}}^{(y)})_{i,j} := (m_{G}^{(y - y_j)})_{i,j}
            \end{array}
        \end{aligned}$$ fulfill the commutative diagram: $$\begin{tikzcd}
            {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{s}}) \arrow[r, "{L}"] \arrow[rd, "\hat{L}"]
            & {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{t}}) \arrow[r, "G"] \arrow[d, "{\mathbb{T}}" ]
            & {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{u}})             \\
            &{
    \mathbcal{L}}(T_{\mathcal{A}}^{\hat{{\mathfrak{t}}}}) \arrow[ur,"\hat{G}"]
        \end{tikzcd}$$

The natural isomorphism between the two corresponding function spaces is given by $$\begin{aligned}
            {\mathbb{T}}: {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{t}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{\hat{{\mathfrak{t}}}}),\  f=(f_1,\ldots,f_m)\mapsto (f_1(\cdot - y_1),\ldots,f_m(\cdot-y_m)) = {\mathbb{T}}f
        \end{aligned}$$ as $f$ and $({\mathbb{T}}f)$ describe the same value distribution on the crystal.[^9] Again, a straightforward calculation yields $$\begin{aligned}
            \begin{array}{rcl}
                [({\mathbb{T}}Lf)(x)]_i & =& [(Lf)(x-y_i)]_i                                                                                              \ =\  \displaystyle\sum\limits_{y \in T_{{\mathcal{A}}}} \displaystyle\sum\limits_{j=1}^{|{\mathfrak{s}}|} (m_L^{(y)})_{i,j}  f_j(x-y_i+y)  \\                 & =&   \displaystyle\sum\limits_{y \in T_{{\mathcal{A}}}} \displaystyle\sum\limits_{j=1}^{|{\mathfrak{s}}|}  (m_L^{(y+y_i)})_{i,j}  f_j(x+y)                 \ =\   [\displaystyle\sum\limits_{y \in T_{{\mathcal{A}}}}  m_{\hat{L}}^{(y)} f(x+y)]_i.
            \end{array}
        \end{aligned}$$ Analogously we find $[( G {\mathbb{T}}^{-1} g)(x)]_i = [\sum_{y \in T_{{\mathcal{A}}}}  (m_{\hat{G}}^{(y)})  g(x+y)]_i$.

Finally, we show that permutations of the entries of the structure element result in a transformation of the non-zero multipliers by corresponding permutation matrices.

\[thm:permuted\_structure\_element\] Consider the two multiplication operators $ L: {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{s}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{t}}})$ and $G: {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{t}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{{\mathfrak{u}}})$ defined by $$\begin{aligned}
            (Lf)(x) = \sum_{y \in T_{{\mathcal{A}}}} m_L^{(y)} f(x+y), \quad m_L^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{t}}|\times |{\mathfrak{s}}|},  \\
            (Gg)(x) = \sum_{y \in T_{{\mathcal{A}}}} m_G^{(y)} g(x+y), \quad m_G^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{u}}| \times |{\mathfrak{t}}|}.
        \end{aligned}$$ Let further $\hat{{\mathfrak{t}}}$ be a structure element which is a permuted version of ${\mathfrak{t}}$, i.e., $$\begin{aligned}
            \hat{{\mathfrak{t}}} & =  (\hat{{\mathfrak{t}}}_1,\ldots,\hat{{\mathfrak{t}}}_m)  =  ({\mathfrak{t}}_{\pi(1)},\ldots, {\mathfrak{t}}_{\pi(m)}) = m_\pi {\mathfrak{t}}\end{aligned}$$ where $m=|{\mathfrak{t}}|$, $\pi:\{1,\ldots,m\} \rightarrow \{1,\ldots,m\}$ is a permutation and $m_\pi \in \{0,1\}^{m\times m}$ the corresponding permutation matrix. Then, the operators $\hat{L}$ and $\hat{G}$ given by $$\begin{aligned}
            \begin{array}{rcllcl}
                (\hat{L}f)(x) &=& \displaystyle\sum\limits_{y \in T_{{\mathcal{A}}}} m_{\hat{L}}^{(y)} f(x+y), & m_{\hat{L}}^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{t}}|\times |{\mathfrak{s}}|} &\text{ with }& m_{\hat{L}}^{(y)} := m_\pi m_{L}^{(y)},\\
                (\hat{G}f)(x) &=& \displaystyle\sum\limits_{y \in T_{{\mathcal{A}}}} m_{\hat{G}}^{(y)} f(x+y), & m_{\hat{G}}^{(y)} \in {
    \mathbb{{\mathcal{C}}}} ^{|{\mathfrak{u}}|\times |{\mathfrak{t}}|} &\text{ with }& m_{\hat{G}}^{(y)} := m_{G}^{(y)} m_\pi^{-1}
            \end{array}
        \end{aligned}$$ fulfill the commutative diagram: $$\begin{tikzcd}
            {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{s}}) \arrow[r, "{L}"] \arrow[rd, "\hat{L}"]
            & {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{t}}) \arrow[r, "G"] \arrow[d, "p" ]
            & {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{u}})             \\
            &{
    \mathbcal{L}}(T_{\mathcal{A}}^{\hat{{\mathfrak{t}}}}) \arrow[ur,"\hat{G}"]
        \end{tikzcd}$$

Due to the fact that the natural isomorphism $p: {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{t}}) \rightarrow  {
    \mathbcal{L}}(T_{\mathcal{A}}^{\hat{{\mathfrak{t}}}})$ is a multiplication operator defined by $(p f)(x) = m_\pi f(x)$ for all $x\in{
    \mathbb{L}}({\mathcal{A}})$, the statement is true due to the rules of computation in \[lem:calc\_rules\_mult\_op\].

and \[thm:lcm\_mat,thm:crystal\_congruence,thm:operator\_sublattice\_coarsening,thm:shifted\_structure\_element,thm:permuted\_structure\_element\] allow for the automatic adjustment of crystal representations within the LFA. The corresponding detailed algorithms which make use of these results are given in \[sec:algorithms\].

Application {#sec:Application}
===========

Before demonstrating the application of aLFA to some selected examples let us briefly recapitulate the individual parts of the framework. The introduction of crystal structures in \[sec:lattices\_crystals\] allow for a canonical description of translationally invariant operators introduced in \[sec:operators\]. Combined with the definition of the dual of a crystal structure in \[def:dual\_lattice\] and a corresponding orthonormal basis in \[thm:wavefunction\_ONB\_vec\], the symbol of any single multiplication operator that is translationally invariant with respect to an arbitrary lattice structure, can be expressed (cf. \[def:crystal\_op\_symbol\]). This combination of choice of basis functions and the explicit connection of operators to their underlying structure enables an automated mixing analysis. This part of the framework can thus be seen on one hand as a unification of positional approaches to LFA by introducing structures that allow for the native treatment of arbitrary translational invariances and on the other hand as a general strategy for the automation of the frequency part back-end of LFA. In order to deal with compositions of operators as encountered in the analysis of iterative methods, the tools provided in \[sec:crystal\_reps\_and\_isos\] allow for an automatic transformation of the underlying crystal structures into compatible representations and the corresponding transformations of the operators. Thus, the only task remaining for the user is to provide any, i.e., the simplest or most convenient, crystal representation of the operators. The following examples show how such a construction can be carried out and as such serve as a tutorial for the use of the algorithms in \[sec:algorithms\]. Annotated Jupyter Notebooks of the presented examples can be found in [@aLFA_NK].

Multicolored block smoother for the tight-binding Hamiltonian of graphene {#subsec:multicolor_graphene}
-------------------------------------------------------------------------

In [@KahlKint2018] a multigrid method for the tight-binding Hamiltonian of the carbon allotrope graphene based on Kaczmarz smoothing is constructed and analyzed using conventional LFA. Due to the hexagonal structure of graphene, the lexicographic ordering of Kaczmarz and the mixing analysis of the coarse grid correction which involved a mixing of eight frequencies, the analysis turned out to be quite lengthy. In this subsection we now want to analyze a two-grid method for this problem where we replace Kaczmarz by an overlapping colored Gauss-Seidel method that allows for better parallelism in the application of the multigrid method. Thus, the goal of this example is two-fold, first we want to show that the tight-binding Hamiltonian can be easily expressed using the native crystal structure of graphene and second that even the analysis of an overlapping block smoother can be carried out with ease using aLFA due to the fact that only a representation of the involved operators is needed.

The graphene structure can be described as a crystal ${
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}})$ where the underlying lattice is triangular, i.e., any three nearby lattice points form an equilateral triangle. We have $$\begin{aligned}
        {
    \mathbb{L}}^{\mathfrak{s}}({\mathcal{A}}),\quad {\mathcal{a}}_1 = (\frac{3}{2}, \frac{\sqrt{3}}{2}) a,\quad {\mathcal{a}}_2 = (\frac{3}{2},-\frac{\sqrt{3}}{2}) a,
    \end{aligned}$$ with the structure element $$\begin{aligned}
        {\mathfrak{s}}= ({\mathfrak{s}}_1,{\mathfrak{s}}_2), \quad {\mathfrak{s}}_1 = ({\mathcal{a}},0)=\frac{1}{3}({\mathcal{a}}_1+{\mathcal{a}}_2), \quad {\mathfrak{s}}_2=(2{\mathcal{a}},0)=\frac{2}{3}({\mathcal{a}}_1+{\mathcal{a}}_2).
    \end{aligned}$$ The parameter $a=  1.42$Å which denotes the distance of two neighboring atoms can be omitted as it does not occur in the tight-binding formulation. An illustration of the crystal has already been given in \[fig:crystal\_example\].

The (nearest neighbor) tight-binding Hamiltonian is a multiplication operator $$\begin{aligned}
        L: {
    \mathbcal{L}}(T_{{\mathcal{A}}} ^{\mathfrak{s}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}}} ^{\mathfrak{s}}), \quad (Lf)(x) := \sum_{y \in T_{{\mathcal{A}}} } m_{L}^{(y)} f(x+y)
    \end{aligned}$$ with non-zero multipliers

\(c) at (0,0) [$=$]{}; ; ;

\(c) at (-,) [$=$]{}; ; ;

\(c) at (,) [$=$]{}; ; ;

\(c) at (-,-) [$=$]{}; ; ;

\(c) at (,-) [$=$]{}; ; ;

as illustrated in \[fig:tbh\_stencil\_and\_4color\_hexagon\], a).

#### Overlapping Hexagons

We now present an overlapping block smoother for the tight-binding Hamiltonian $L$. Consider the non-disjoint splitting (or coloring) of the crystal into hexagons as depicted in \[fig:tbh\_stencil\_and\_4color\_hexagon\], b). This splitting has a translational invariance of ${\mathcal{C}}= 2{\mathcal{A}}$. Rewriting the graphene crystal $T_{\mathcal{A}}^{\mathfrak{s}}$ with respect to this coarser lattice ${
    \mathbb{L}}({\mathcal{C}})$, we find $T_{\mathcal{A}}^{\mathfrak{s}}\cong T_{\mathcal{C}}^{\mathfrak{t}}$ with the structure element $$\begin{aligned}
        {\mathfrak{t}}=({\mathfrak{t}}_1,\ldots,{\mathfrak{t}}_8)= ({\mathfrak{s}}, {\mathfrak{s}}+{\mathcal{a}}_1, {\mathfrak{s}}+{\mathcal{a}}_2, {\mathfrak{s}}+{\mathcal{a}}_1+{\mathcal{a}}_2).
    \end{aligned}$$ The splitting is then given by the structure elements $$\begin{aligned}
        {\mathfrak{t}}^{(1)}&= ({\mathfrak{t}}_1^{(1)},\ldots,{\mathfrak{t}}_6^{(1)})  = ({\mathfrak{t}}_2,{\mathfrak{t}}_3,{\mathfrak{t}}_4,{\mathfrak{t}}_5,{\mathfrak{t}}_6,{\mathfrak{t}}_7), \\
        {\mathfrak{t}}^{(2)}&= {\mathfrak{t}}^{(1)} + {\mathcal{a}}_1,\         {\mathfrak{t}}^{(3)}=  {\mathfrak{t}}^{(1)} + {\mathcal{a}}_2,\text{ and }         {\mathfrak{t}}^{(4)}=  {\mathfrak{t}}^{(1)} + {\mathcal{a}}_1 + {\mathcal{a}}_2,     \end{aligned}$$ such that $T_{\mathcal{C}}^{{\mathfrak{t}}^{(1)}} \widehat{=} \CGcaptionAab[]$, $T_{\mathcal{C}}^{{\mathfrak{t}}^{(2)}} \widehat{=} \CGcaptionAbb[]$, $T_{\mathcal{C}}^{{\mathfrak{t}}^{(3)}} \widehat{=} \CGcaptionAba[]$ and $T_{\mathcal{C}}^{{\mathfrak{t}}^{(4)}} \widehat{=} \CGcaptionAaa[]$. In the case of the nearest neighbor Hamiltonian any two hexagons of the same color do not interact. Thus, a relaxed sweep on the unknowns of one color is cheap and, in addition, yields a good degree of parallelism.

There are two options to describe the multiplication operators of the multicolored block smoother. On the one hand, we can construct them from scratch by defining the multiplication operator structure as well as the non-zero multipliers. On the other hand, we can derive the structure and non-zero multipliers of the operators corresponding to the colored blocking from the tight-binding Hamiltonian $L$ by exploiting the fact that a colored block corresponds to the operator $L$ restricted to this block. In order to proceed with the latter approach, we first need to find the description $\hat{L}$ of the underlying operator, the tight-binding Hamiltonian $L$, with respect to the translational invariance of the splitting ${\mathcal{C}}$. After that the structure element needs to be adjusted such that all unknowns and their couplings among each other are found within the central multiplier $m_{\hat{L}}^{(0)}$. Consider the structure element ${\mathfrak{t}}$. As can be seen in \[fig:tbh\_stencil\_and\_4color\_hexagon\], b), the coupling among the unknowns ${\mathfrak{t}}^{(\ell)}$ which we want to update simultaneously are found in the multipliers:

[c|c|c|c]{} ${\mathfrak{t}}^{(1)}$ & ${\mathfrak{t}}^{(2)}$ & ${\mathfrak{t}}^{(3)}$ & ${\mathfrak{t}}^{(4)}$\

------------------------------------------------------------------------

$m_{\hat{L}}^{(0)}$ & $m_{\hat{L}}^{(0)},\ m_{\hat{L}}^{(\pm {\mathcal{c}}_1)}$ & $m_{\hat{L}}^{(0)},\ m_{\hat{L}}^{(\pm{\mathcal{c}}_2)}$ & $m_{\hat{L}}^{(0)},\ m_{\hat{L}}^{(\pm{\mathcal{c}}_1)},\ m_{\hat{L}}^{(\pm{\mathcal{c}}_2)},\  m_{\hat{L}}^{(\pm({\mathcal{c}}_1-{\mathcal{c}}_2))}$  

Thus, in order to obtain suitable descriptions for ${\mathfrak{t}}^{(\ell)}$, $\ell\in\{2,3,4\}$, we need to consider shifted versions ${\mathfrak{t}}+ \tau^{(\ell)}$. For example $$\label{eq:shiftedhex}
        \begin{bmatrix} \tau^{(1)}, & \tau^{(2)}, & \tau^{(3)}, & \tau^{(4)}\end{bmatrix} := \begin{bmatrix} 0, & {\mathcal{a}}_1, & {\mathcal{a}}_2, & {\mathcal{a}}_1 + {\mathcal{a}}_2 \end{bmatrix} .$$ Finally, the error propagator corresponding to a single color relaxed block Gauss-Seidel update can be written as $$\begin{aligned}
        G^{(\ell)} = (I - \omega (S^{(\ell)})^\dagger \hat{L}) : {
    \mathbcal{L}}(T_{\mathcal{C}}^{{\mathfrak{t}}+ \tau^{(\ell)} }) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{C}}^{{\mathfrak{t}}+ \tau^{(\ell)} })
    \end{aligned}$$ with $(S^{(\ell)}f)(x) := (m_P m_{\hat{L}}^{(0)} m_P)f(x)$ where $m_P$ is the diagonal matrix $$\label{eq:mpii_hex}
        (m_P)_{ii} = \begin{cases} 1 & \text{ for }i\in\{2,3,4,5,6,7\} \\
            0 & \text{ else}.
        \end{cases}$$

We summarize the algorithmical steps to obtain the error propagators. The error propagator of a successive update of the four colors corresponds to the product of the error propagators $G^{(\ell)} = ( I - \omega (S^{(\ell)})^\dagger \hat{L}^{(\ell)})$. Now that all operators of the overlapping colored block Gauss-Seidel smoother are defined, we can compute its spectrum. The computation of the eigenvalues of the error propagator is carried out with \[alg:compute\_spectral\_radii\] via $$\textproc{ComputeSpectrum}(g,I,L,S^{(1)},\ldots,S^{(4)}).$$ The function $g$ denotes the composition of the error propagators $$\begin{aligned}
(I,L,S^{(1)},\ldots,S^{(4)}) \xmapsto{\enskip g\enskip } \prod_{\ell=1}^4 ( I - \omega (S^{(\ell)})^\dagger L)=:G.
    \end{aligned}$$ Note, that in this algorithm all operators are checked for compatibility and any incompatibility with respect to domains and codomains is dealt with using the transformations introduced in \[sec:crystal\_reps\_and\_isos\]. For $\omega=\frac{1}{2}$ we obtain the plot in \[fig:hexagon\_smooth\_twogrid\], a). As the largest spectral radius is greater than one, this method cannot be used as a standalone solver.

#### Coarse grid correction

In [@KahlKint2018] a Galerkin coarse grid correction is used with a corresponding coarse grid correction operator that is defined by $$\begin{aligned}
        E = (I- P(L_c)^\dagger R L)
    \end{aligned}$$ with $L_c = RLP$ and $P=R^T$. Just as in the description of the smoother, the shift invariance of the coarse grid is ${\mathcal{C}}=2{\mathcal{A}}$. Defining the coarse structure element by $2{\mathfrak{s}}$, and with ${\mathfrak{t}}$ denoting the structure element of the fine crystal according to the coarse lattice basis, the restriction operator can be described as $$\begin{aligned}
        R : {
    \mathbcal{L}}(T^{\mathfrak{t}}_{2{\mathcal{A}}}) \rightarrow {
    \mathbcal{L}}(T^{2{\mathfrak{s}}}_{2{\mathcal{A}}}).
    \end{aligned}$$ The multipliers of the restriction operator according to [@KahlKint2018] then read

This coarse grid correction is constructed in a way, such that it preserves the kernel of the tight-binding Hamiltonian $L$, where with ${\mathcal{B}}=\mat{{\mathcal{b}}_1 & {\mathcal{b}}_2} = {\mathcal{A}}^{-T}$, $$\operatorname{ker}(L)=\operatorname{span}{\{\mat{{e^{2\pi i{\langle K,x\rangle_{2}}}} \\ 0}, \mat{0 \\ {e^{2\pi i{\langle K,x\rangle_{2}}}}} \, : \, K \in\{\frac{1}{3} {\mathcal{b}}_1 + \frac{2}{3} {\mathcal{b}}_2,\frac{2}{3} {\mathcal{b}}_1 + \frac{1}{3} {\mathcal{b}}_2\} \}}.$$ The frequencies corresponding to the kernel modes are known as the *Dirac points*. The two-grid analysis can now be carried out using \[alg:compute\_spectral\_radii\] via $$\text{\textproc{ComputeSpectrum}}(f,I,L,S^{(1)},\ldots,S^{(4)}, R).$$ The function $f$ denotes the composition of the two-grid error propagator $$\begin{aligned}
(I,L,S^{(1)},\ldots,S^{(4)},R) \xmapsto{\enskip f\enskip } GEG.    \end{aligned}$$

A plot of the spectral radii of the two-grid error propagator is given in b) of \[fig:hexagon\_smooth\_twogrid\] which shows that the two-grid method converges with a convergence rate of $\rho_{\max} \approx 0.167$. Thus, this new method with overlapping colored block Gauss-Seidel smoothing not only yields opportunities for parallel computations, but also converges faster than the old approach which used Kaczmarz smoothing.

In order to double-check the results of the developed theory, we show in \[tbl:graphene\_twogrid\_results\] that the asymptotic convergence rate of the two-grid method with random initial guess $x_0$ and right-hand-side $0$ coincides with the convergence rate obtained in the LFA with a relative accuracy of roughly $.002\%$. This comes as no surprise as the theory is exact for this problem with periodic boundary conditions and we have chosen the sampling of the frequency space in accordance with the problem size (cf. \[rem:freq\_sampling\]).

   iteration $i$   $||r_i||_2:=||b - A x_i||_2$   $\rho_i:=\frac{||r_i||_2}{||r_{i-1}||_2}$   $\frac{|\rho_i - \rho_{\text{analytic}}|}{\rho_{\text{analytic}}}$
  --------------- ------------------------------ ------------------------------------------- --------------------------------------------------------------------
        398               1.161369e-312                          0.16685534                                                0.00220%
        399               1.937806e-313                          0.16685539                                                0.00217%
        400               3.233335e-314                          0.16685545                                                0.00213%

  : Convergence history of the two-grid method applied to the tight-binding Hamiltonian of graphene with $(41\times 41)\cdot 2^3$ unknowns/atoms and periodic boundary conditions. The reported asymptotic convergence rate $\rho_i$ coincides with high precision to the convergence estimate $\rho_{\text{analytic}} = \sup_{k \in T^*_{2{\mathcal{A}},41\cdot 2{\mathcal{A}}}} {\{| \lambda | \, : \, \lambda \text{ eigenvalue of } G_kE_kG_k\}} = 0.16685901$ obtained in aLFA.[]{data-label="tbl:graphene_twogrid_results"}

Two-level analysis for the curl-curl equation {#subsec:curlcurl}
---------------------------------------------

In [@BoonLentVand2008] a complete two-level analysis is carried out for the $2$-dimensional curl-curl formulation of Maxwell’s equations using conventional LFA. In this subsection we want to reproduce the results of the two-grid method discussed in this paper making use of the native crystal structure of the staggered discretization of the curl-curl equations. The method consists of a so-called half-hybrid smoother, introduced in [@Hipt1998], and a Galerkin coarse grid correction. The degrees of freedom of the discrete curl-curl equation $$\begin{aligned}
\label{eq:discrete_curl}
        (\frac{1}{h^2}K_{cc} + \sigma M)x= b
    \end{aligned}$$ are associated with the edges of a quadrilateral lattice ${
    \mathbb{L}}(\frac{1}{h} {\mathcal{A}})$ with $${\mathcal{A}}= \mat{{\mathcal{a}}_1 &  {\mathcal{a}}_2}= \mat{1 & 0 \\ 0 & 1}.$$ By multiplying \[eq:discrete\_curl\] with $h^2$, the grid size can be incorporated into the material parameter $\sigma_h := h^2\sigma$. Then the discrete operator $K=K_{cc} + \sigma_h M$ can be expressed as a multiplication operator by $$\begin{aligned}
        K: {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{e}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{e}}) \quad \text{with} \quad 
        {\mathfrak{e}}= ({\mathfrak{e}}_h,{\mathfrak{e}}_v) = (\frac{1}{2}{\mathcal{a}}_1, \frac{1}{2}{\mathcal{a}}_2),
    \end{aligned}$$ where the elements ${
    \mathbb{L}}({\mathcal{A}}) + {\mathfrak{e}}_h$ and ${
    \mathbb{L}}({\mathcal{A}}) + {\mathfrak{e}}_v$ correspond to the (midpoints of the) horizontal and the vertical edge, respectively (cf. \[fig:curlcurl\_crystal\]).

According to [@BoonLentVand2008] the multipliers are given by:

#### Hybrid smoother

The (half) hybrid smoother uses an auxiliary space correction. That is, it can be seen as a two-grid method itself where smoothing replaces the exact inversion of the auxiliary space operator. The construction of this operator is as follows.

As already stated, we want to reproduce the results in [@BoonLentVand2008]. Due to the fact that lexicographic Gauss-Seidel smoothers are used, we introduce the exact same ordering of the lattice points used in this paper, that is, the coordinates $y\in {
    \mathbb{L}}({\mathcal{A}})$ are ordered from bottom to top and left to right which leads us to the following definition for $x,y\in{
    \mathbb{L}}({\mathcal{A}})$: $$x=(x_1,x_2) < y=(y_1,y_2) :\Leftrightarrow x_2 < y_2 \text{ or } \left(x_2 = y_2 \text{ and } x_1 < y_1\right) .$$

The error propagator of the smoother on the original crystal $T_{\mathcal{A}}^{\mathfrak{e}}$ is a node-based lexicographic Gauss-Seidel iteration. When updating a horizontal edge $x + {\mathfrak{e}}_h$, $x\in{
    \mathbb{L}}({\mathcal{A}})$ it is assumed that all edges $y+{\mathfrak{e}}_h$ and $y+\hat{{\mathfrak{e}}}_v$, $\hat{{\mathfrak{e}}}_v:={\mathfrak{e}}_v+{\mathcal{a}}_1 - {\mathcal{a}}_2$, with $y<x$ are already updated. When updating the edge $y+\hat{{\mathfrak{e}}}_v$, it is additionally assumed that $x+{\mathfrak{e}}_h$ is already updated (cf. [@BoonLentVand2008 Figure 6.1]). In order to obtain the error propagator which exactly represents this ordering, it is convenient to rewrite the operator $K$ with respect to this representation of the structure element, that is $$\begin{aligned}
            K \cong \hat{K}: {
    \mathbcal{L}}(T_{\mathcal{A}}^{\hat{{\mathfrak{e}}}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{\hat{{\mathfrak{e}}}}) \quad \text{with} \quad 
            \hat{{\mathfrak{e}}} = ({\mathfrak{e}}_h,\hat{{\mathfrak{e}}}_v) \cong {\mathfrak{e}},
    \end{aligned}$$ which can be obtained with \[alg:struct\_elem\_change\] via $$\begin{aligned}
            \hat{K}=\text{\textproc{ChangeStructureElement}}(K,\hat{{\mathfrak{e}}},\hat{{\mathfrak{e}}}).        \end{aligned}$$ The corresponding error propagator is then given by $$\begin{aligned}
        G_E = (I - S_E^\dagger K ): {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{e}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{e}})
    \end{aligned}$$ with the multipliers $$\begin{aligned}
        m_{S_E}^{(y)} = 
        \begin{cases}
            m_{\hat{K}}^{(y)} & \text{if } y < 0, \\
            \texttt{tril}(m_{\hat{K}}^{(y)}) & \text{if } y=0, \\
            0 & \text {else.}
        \end{cases}
    \end{aligned}$$ In here, only the lower triangular part $\texttt{tril}(m_{\hat{K}}^{(0)})$ of the central multiplier is used due to the fact that it represents a Gauss-Seidel sweep where a horizontal edge is updated before a vertical edge.

The crystal, where the auxiliary space system is formulated, is given by $T_{\mathcal{A}}^{(0)}$, i.e., the nodal points between the edges. The transfer operator to this crystal is the discrete gradient operator defined by $R_N: {
    \mathbcal{L}}(T_{{\mathcal{A}}}^{\mathfrak{e}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}}}^{(0)})$ with non-zero multipliers $$\begin{aligned}
        \begin{array}{rccrcc}
            m_{R_N}^{(-{\mathcal{a}}_1)} &=& \mat{1 & 0} & m_{R_N}^{(0)} &=& \mat{-1 & -1} \\[1em]
            & & & m_{R_N}^{(-{\mathcal{a}}_2)} &= &\mat{0 & 1}.
        \end{array}
    \end{aligned}$$ Using a Galerkin construction, i.e., $P_N = R_N^T$, the coarse grid operator $K_N = R_N K P_N$ can be obtained using the computation rules in \[lem:calc\_rules\_mult\_op\]. Inversion of this auxiliary space operator is then approximated by Gauss-Seidel. Thus, the auxiliary space correction with a single smoothing step on the auxiliary space is given by $$\begin{aligned}
        G_N = (I - P_N S_N^\dagger R_N K) : {
    \mathbcal{L}}( T_{\mathcal{A}}^{\mathfrak{e}}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{\mathfrak{e}}),\end{aligned}$$ where $S_N: {
    \mathbcal{L}}(T_{\mathcal{A}}^{(0)}) \rightarrow {
    \mathbcal{L}}(T_{\mathcal{A}}^{(0)})$ consists of the scalar multipliers $$\begin{aligned}
        m_{S_N}^{(y)} = \begin{cases} m_{K_N}^{(y)} & \text{if } y \leq 0, \\
            0 & \text {else.}\end{cases}
    \end{aligned}$$ With this definition of the auxiliary space correction, the half-hybrid smoother consists of the following steps.

1.  Smooth on $Kx=b$: $x \leftarrow (I+S_E^\dagger K)x + S_E^\dagger b$,

2.  Restrict the residual: $r_N \leftarrow R_N(b-Kx)$,

3.  Smooth on $K_N x_N= r_N$ with zero initial guess: $x_N \leftarrow S_N^\dagger r_N$,

4.  Prolongate to the primal crystal $x \leftarrow x + P_N x_N$.

The smoother is analyzed analogously to \[subsec:multicolor\_graphene\] with \[alg:compute\_spectral\_radii\] by computing the spectral radii via $$\text{\textproc{ComputeSpectrum}}(g,I,K,S_E,R_N,S_N),$$ where $g$ denotes the composition of the error propagators $$(I,K,S_E,R_N,S_N) \xmapsto{\enskip g\enskip } G_N G_E =:G.$$ In \[fig:LFA\_curlcurl\] a) a contour plot of $\rho(G_k)$ with $\sigma_h = 0.01$, with respect to $k\in {\mathcal{A}}^{-T}[ -\frac{1}{4},\frac{3}{4})$ is given. It corresponds to the result given in [@BoonLentVand2008 Figure 6.2, right].

#### Coarse grid correction

In [@BoonLentVand2008] a Galerkin coarse grid correction is used corresponding to the error propagator $$\begin{aligned}
        E = (I- P(K_c)^\dagger R K)
    \end{aligned}$$ with $K_c = RKP$ and $P=R^T$. In here, the coarse crystal is ${
    \mathbb{L}}^{2{\mathfrak{e}}}(2{\mathcal{A}})$ and the original crystal ${
    \mathbb{L}}^{{\mathfrak{e}}}({\mathcal{A}})$ with respect to the lattice $2{\mathcal{A}}$ is given by ${
    \mathbb{L}}^{{\mathfrak{f}}}(2{\mathcal{A}})$ with $$\begin{aligned}
        {\mathfrak{f}}= ({\mathfrak{e}}, {\mathfrak{e}}+ {\mathcal{a}}_1, {\mathfrak{e}}+ {\mathcal{a}}_2, {\mathfrak{e}}+ {\mathcal{a}}_1 + {\mathcal{a}}_2).
    \end{aligned}$$ Then, according to [@BoonLentVand2008] the restriction operator is given as $R : {
    \mathbcal{L}}(T^{\mathfrak{f}}_{2{\mathcal{A}}}) \rightarrow {
    \mathbcal{L}}(T^{2{\mathfrak{s}}}_{2{\mathcal{A}}})$ with the multipliers:

Here, each coarse horizontal and vertical edge is connected to its six nearest edges of the same type (cf. \[fig:curlcurl\_crystal\] b)). The prolongation and coarse grid operator $P$ and $K_c$ can be obtained via \[lem:calc\_rules\_mult\_op\].

The spectral radii of the two-grid method can then be obtained via $$\begin{aligned}
                        \text{\textproc{ComputeSpectrum}}(f,I,K,S_E,R_N,S_N,R)             \end{aligned}$$ where $f$ denotes the composition of the two-grid error propagator $$\begin{aligned}
                        (I,K,S_E,R_N,S_N,R) &\xmapsto{\enskip f\enskip } G E G =: M.
                    \end{aligned}$$ b) shows the spectral radii, $\sup_k \rho(M_k)$, of the two-grid error propagator as a function of $\sigma_h$. We used exactly the same orderings in pre- and post-smoothing, as described in [@BoonLentVand2008]. Furthermore, we double-checked the results with convergence tests similar to \[tbl:graphene\_twogrid\_results\]. However, the obtained results show major differences to the one in [@BoonLentVand2008 Figure 8.1, left], which leads us to believe that our assumptions on the lexicographic orderings employed in pre- or post-smoothing in [@BoonLentVand2008] are wrong, as these have a large impact on the convergence rates.

Conclusion
==========

In this paper we present an LFA framework that is applicable to arbitrary crystal structures, which are encountered in many applications, e.g., systems of partial differential equations, block smoothers or tight-binding formulations. This was achieved by introducing a rigorous notion of crystal structures and translationally invariant operators that manipulate value distributions on crystal structures. Based on these structures we introduced a complete framework to modify and transform these operators with respect to different crystal representations by using normal forms of integer linear algebra. This allowed us to automate the LFA in both position space, i.e., in terms of multiplication operators/stencils, and frequency space, i.e., in terms of canonical basis functions and samplings. In two examples we showed that the approach can be used for complicated operators, i.e., hexagonal grids, overlapping block smoothers and hybrid smoothers, without requiring insight into the frequency space back-end of LFA. An explicit mixing calculation is no longer needed, and the user only has to provide any representation of the individual operators. The transformation of the individual operators to a compatible representation and the subsequent frequency analysis is then carried out automatically. Even though we have limited ourselves in the examples in this paper to $2$-dimensional problems our automated LFA can be applied to operators in higher dimension, where the difference is a larger set of primitive vectors defining the underlying lattice.

The automation presented in this paper does have some limitation. Each individual operator in the analysis is only allowed to change each value of the value distribution at most once. This limitation solely restricts the class of smoothers that can be analyzed with this approach. Any sequential, i.e., lexicographic, smoother with overlapping update regions changes values in the overlap multiple times in one application. This cannot be easily translated to a corresponding local multiplication operator, but it can be dealt with in frequency space (cf. [@MacLOost2011; @Molenaar1991; @Sivaloganathan:1991:ULM]). This particular treatment of sequential overlap is momentarily not covered in our framework. Note, that the mere presence of overlap is not the problem here. By introducing a coloring, such that the complete sweep can be split into a sequence of updates where each one of them only changes values at most once, automated LFA can be applied (cf. \[subsec:multicolor\_graphene\]). Due to the fact that a coloring in overlapping approaches also favors parallelism over their sequential counterparts, we feel that this limitation is relatively minor when targeting actual applications. An open-source implementation of the automated LFA framework [@aLFA_NK] is freely available on GitLab[^10]. Jupyter Notebooks for all examples from \[sec:example,subsec:multicolor\_graphene,subsec:curlcurl\] are included in this software package as well.

Rules of computation {#sec:rules_of_comp}
====================

Calculus of multiplication operators plays a key-role in local Fourier analysis. In this section we list all elementary operations, such as addition and multiplication. Proofs for these rules can be obtained by straightforward calculation.

\[lem:calc\_rules\_mult\_op\] Let two multiplication operators be given by $$\begin{aligned}
            \begin{array}{rcl}
                L:{
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{s}}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{t}}}), &
                (Lf)(x) = \displaystyle\sum_{y\in T_{{\mathcal{A}},{\mathcal{Z}}}}  m_L^{(y)} f(x+y),& m_L^{(y)} \in{
    \mathbb{C}}^{|{{\mathfrak{t}}}|\times |{{\mathfrak{s}}}|},\\
                G:{
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{u}}})\rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{v}}}),&                 (Gf)(x) = \displaystyle\sum_{y\in T_{{\mathcal{A}},{\mathcal{Z}}}}  m_G^{(y)} f(x+y), & m_G^{(y)} \in{
    \mathbb{C}}^{|{{\mathfrak{v}}}|\times |{{\mathfrak{u}}}|}.
            \end{array}
        \end{aligned}$$ Then the following operators are multiplication operators as well:

1.  If ${\mathfrak{s}}= {\mathfrak{u}}$ and ${\mathfrak{t}}={\mathfrak{v}}$, then $
                L+G:{
        \mathbcal{L}}(T_{{
        \mathcal{A}},{
        \mathcal{Z}}}^{\mathfrak{s}}) \rightarrow {
        \mathbcal{L}}(T_{{
        \mathcal{A}},{
        \mathcal{Z}}}^{\mathfrak{t}}) \text{ with } m_{L+G}^{(y)} =  m_L^{(y)}+m_G^{(y)}.$

2.  If ${\mathfrak{v}}= {\mathfrak{s}}$, then $LG:{
        \mathbcal{L}}(T_{{
        \mathcal{A}},{
        \mathcal{Z}}}^{{\mathfrak{u}}}) \rightarrow {
        \mathbcal{L}}(T_{{
        \mathcal{A}},{
        \mathcal{Z}}}^{{{\mathfrak{t}}}}) \text{ with } m_{LG}^{(z)} = \sum\limits_{y+w=z} m_L^{(y)} \cdot m_G^{(w)}.$

3.  The adjoint is given by $L^*:{
        \mathbcal{L}}(T_{{
        \mathcal{A}},{
        \mathcal{Z}}}^{\mathfrak{t}}) \rightarrow {
        \mathbcal{L}}(T_{{
        \mathcal{A}},{
        \mathcal{Z}}}^{\mathfrak{s}}) \text{ with }m_{L^*}^{(y)} = (m_{L}^{(-y)})^*.$

\[thm:Eigenvalues\_Crystal\_Mult\_Op\] Let two multiplication operators be given by $$\begin{aligned}
            \begin{array}{rcl}
                L:{
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{s}}}) \rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{t}}}), &
                (Lf)(x) = \displaystyle\sum_{y\in T_{{\mathcal{A}},{\mathcal{Z}}}}  m_L^{(y)} f(x+y),& m_L^{(y)} \in{
    \mathbb{C}}^{|{{\mathfrak{t}}}|\times |{{\mathfrak{s}}}|},\\
                G:{
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{u}}})\rightarrow {
    \mathbcal{L}}(T_{{\mathcal{A}},{\mathcal{Z}}}^{{\mathfrak{v}}}),&                 (Gf)(x) = \displaystyle\sum_{y\in T_{{\mathcal{A}},{\mathcal{Z}}}}  m_G^{(y)} f(x+y), & m_G^{(y)} \in{
    \mathbb{C}}^{|{{\mathfrak{v}}}|\times |{{\mathfrak{u}}}|}
            \end{array}
        \end{aligned}$$ with corresponding symbols $L_k$ and $G_k$. Then we have the following statements.

1.  Assuming that ${{\mathfrak{s}}} = {{\mathfrak{u}}}$ and ${{\mathfrak{t}}} = {{\mathfrak{v}}}$, the symbols of $L+G$ are given by $L_k+ G_k$.

2.  Assuming that ${{\mathfrak{v}}} = {{\mathfrak{s}}}$, the symbols of $L\cdot G$ are given by $L_k\cdot G_k$.

3.  The symbols of $L^*$ are given by $L_k^*$.

4.  The symbols $(L^\dagger)_k$ are given by $(L_k)^\dagger$.

Algorithms {#sec:algorithms}
==========

[ ]{}

[ ]{} ($\triangleright$ See \[def:crystal\_op\_and\_normalform\])[$G = \func{Normalize}{L}$]{}[ ${\mathfrak{u}}_j = {\mathfrak{s}}_j - {\mathcal{A}}\lfloor {\mathcal{A}}^{-1}{\mathfrak{s}}_j \rfloor $ for all $j=1,\ldots,|{\mathfrak{s}}|$ ${\mathfrak{v}}_j = {\mathfrak{t}}_j - {\mathcal{A}}\lfloor {\mathcal{A}}^{-1}{\mathfrak{t}}_j \rfloor $ for all $j=1,\ldots,|{\mathfrak{t}}|$ Sort ${\mathfrak{u}}$ and ${\mathfrak{v}}$ lexicographically $G = \func{ChangeStructureElement}{L,{\mathfrak{u}},{\mathfrak{v}}}$ ]{}

[ ]{} ($\triangleright$ See \[cor:quotient\_latticepointlist\])[${\mathfrak{s}}= \func{ElementsInQuotientSpace}{{\mathcal{A}},{\mathcal{C}}}$]{}[ $H = $ Hermite normal form of ${\mathcal{A}}^{-1}{\mathcal{C}}$ $m = \prod_{i=1}^n h_{i,i}$ ${\mathfrak{s}}_i = 0$ for all $i=0,\ldots,m$ ]{}

[ ]{} ($\triangleright$ See \[thm:operator\_sublattice\_coarsening\])[$G = \func{LatticeCoarsening}{L,{\mathcal{C}}}$]{}[ ${\mathfrak{e}}= \func{ElementsInQuotientSpace}{{\mathcal{A}},{\mathcal{C}}}$ ${\mathfrak{u}}= ({\mathfrak{e}}_1 + {\mathfrak{s}},\ldots, {\mathfrak{e}}_{|{\mathfrak{e}}|} + {\mathfrak{s}}), {\mathfrak{v}}= ({\mathfrak{e}}_1 + {\mathfrak{t}},\ldots, {\mathfrak{e}}_{|{\mathfrak{e}}|} + {\mathfrak{t}})$ $(m_{{G}}^{({y})})=0 \in {
    \mathbb{C}}^{|{\mathfrak{v}}|\times |{\mathfrak{u}}|}$ for all ${y}\in{
    \mathbb{L}}({\mathcal{C}})$ ]{}

[ ]{}

[ ]{} ($\triangleright$ See \[thm:shifted\_structure\_element,thm:permuted\_structure\_element\])[$\hat{L} = \func{ChangeStructureElement}{L,{{\mathfrak{u}}},{{\mathfrak{v}}}}$]{}[ $m_\pi =0 \in \{0,1\}^{|{\mathfrak{s}}|\times|{\mathfrak{s}}| }$, $m_\sigma = 0 \in \{0,1\}^{|{\mathfrak{t}}|\times|{\mathfrak{t}}| }$ ($\triangleright$ Compute changes in ${\mathfrak{d}}$)[$(i,j) \in \{1,\ldots,|{\mathfrak{s}}|\}^2$]{}[ ]{} ($\triangleright$ Compute changes in ${\mathfrak{c}}$)[$(i,j) \in \{1,\ldots,|{\mathfrak{t}}|\}^2$]{}[ ]{} $(m_{\hat{L}}^{({y})})=0 \in {
    \mathbb{C}}^{|{\mathfrak{t}}|\times |{\mathfrak{s}}|}$ for all ${y}\in{
    \mathbb{L}}({\mathcal{A}})$ $m_{\hat{L}}^{(y)} =  m_\sigma \cdot  m_{\hat{L}}^{(y)}\cdot m_\pi^{-1}$ for all ${y}\in{
    \mathbb{L}}({\mathcal{A}})$ ]{}

[ ]{} ($\triangleright$ See \[thm:lcm\_mat\])[${\mathcal{C}}= \func{LeastCommonMultiple}{{\mathcal{A}},{\mathcal{B}}}$]{}[ Find an integer $r$, s.t. $M=r{\mathcal{A}}^{-1}{\mathcal{B}}$ is integral Compute Smith normal form $S=V^{-1}MT^{-1}$ of $M$ $(N_{{\mathcal{B}}})_{i,i} = r \cdot \operatorname{gcd}(r,s_{i})^{-1}$, ${\mathcal{C}}= {\mathcal{B}}T^{-1} N_{{\mathcal{B}}}$ ]{}

[^1]: School of Mathematics and Natural Sciences, University of Wuppertal, 42097 Germany, `{kkahl,kintscher}@math.uni-wuppertal.de`

[^2]: This work was partially funded by Deutsche Forschungsgemeinschaft (DFG) Transregional Collaborative Research Centre 55 (SFB/TRR55)

[^3]: [`github.com/hrittich/lfa-lab`](https://github.com/hrittich/lfa-lab)

[^4]: [`gitlab.com/NilsKintscher/alfa`](https://gitlab.com/NilsKintscher/alfa)

[^5]: All proofs of \[lem:sublattice,thm:integer\_normalforms,thm:lattice\_equiv\_unimod\] can be found in [@Schrijver:1986:TLI:17634].

[^6]: Bounded by the number of lattice points on the (arbitrarily large) torus.

[^7]: Recall that a unique list of representatives $T_{{\mathcal{A}},{\mathcal{C}}} = \{{\mathfrak{t}}_1,\ldots,{\mathfrak{t}}_p\}$ can be found via \[cor:quotient\_latticepointlist\].

[^8]: In case ${\mathfrak{d}}_i={\mathfrak{d}}_j$ or ${\mathfrak{c}}_i = {\mathfrak{c}}_j$ for any $i\neq j$ a consistent ordering of $i,j$ has to be defined a priori.

[^9]: The left-hand side of $[({\mathbb{T}}f)(x)]_i=f_i(x-y_i)$ corresponds to the value at position $x + {\hat{{\mathfrak{t}}}}_i = (x-y_i)+{\mathfrak{t}}_i$ which coincides with the position of the value of the right-hand side.

[^10]: [`gitlab.com/NilsKintscher/alfa`](https://gitlab.com/NilsKintscher/alfa)
