Automated local Fourier analysis (aLFA)

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.


Introduction
Local Fourier analysis (LFA) is a powerful tool used in the construction and analysis of multigrid methods introduced in [4].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 [18].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 [5,15,22].In other cases a similar convergence rate can usually be obtained by additional processing [9,23].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 [23,24].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) [23].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.[23,24]).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 [6] and on hexagonal tilings in [26].Further, it has been applied to edge-based quadriliteral discretizations [3], regular Voronoi meshes associated with acute triangular grids [16], edge-based discretizations on triangular grids [17] and jumping coefficients on rectangular grids [2].These papers usually do a complete two-grid analysis.Three-grid LFA for rectangular grids can for example be found in [25].
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 to a certain degree [14,24].In contrast to the software described in [24], which basically consists of a collection of templates corresponding to certain smoother and restriction/prolongation strategies for specific problems, the software resulting from [14] can be used to describe and analyze arbitrary translationally invariant operators on rectangular grids.Even though this software automates several technical steps of the LFA, there are still tasks which need be done manually.For example in the analysis of block smoothers one has to define the diagonal blocks of the system matrix which are inverted.In general, all operators of the error propagator need to be supplied on a single shared translational invariance of rectangular shape which has to be determined a priori by the user.Furthermore, how to proceed if the blocks of a block smoother do not fit into one rectangular domain, which for example is the case for some overlapping block smoothers, is not discussed in [14].
In this paper we present an approach to LFA which requires minimal user input and fully automates the process of LFA by additional (automated) preprocessing in position space, which in turn lets us fully automate the frequency space back-end as well.For example, in order to analyze a block method, the presented framework requires the user to only determine the unknowns corresponding to each block.In contrast to the approach developed in this paper, the bulk of the analysis in the previously mentioned references is carried out in frequency space.LFA with an emphasize on position space is very rare in the literature.In [8] 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 intro-duce a mathematical framework for the analysis of translationally invariant operators which alter value distributions on lattices and crystals [1,19].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 tightbinding descriptions of solid-state physics [1] 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 nonmatching) descriptions of the underlying repetitive structures, i.e., each operator can be supplied in the simplest or most convenient representation.In contrast to previously developed LFA we explicitly include the connection of the operator to the underlying structure.This allows us on one hand to easily manipulate 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 a full 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 [19].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 (e.g., incorporating overlapping block smoothers) on complicated structures.
The paper is organized as follows.In section 2 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 section 3.After that, section 4 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 section 3. Section 5 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 section 2.1.Using these results we obtain the algorithms given in appendix B which, in combination with the arithmetic of multiplication operators given in appendix A, realize the automated LFA.Finally, section 6 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 [3] as a further illustration of how our approach is applied to complex error propagators and to double check its results.

Lattices and Crystals
In order to be able to automate and generalize the process of LFA we first have to review some basic definitions of integer linear algebra and crystallography [1,19].An (ideal) crystal is an infinite repetition, defined by a lattice, of a structure element.Definition 2.1.Let a 1 , a 2 , . . ., a n ∈ R n be linearly independent.An n-dimensional lattice L is the set of points The vectors a 1 , a 2 , . . ., a n are known as the primitive vectors or lattice basis.Using matrix notation, i.e., A := a 1 a 2 . . .a n , we can abbreviate the notation by Without loss of generality we restrict the definition of the second component of a crystal, the structure element, to primitive cells of the lattice.
volume of space that, if translated by all vectors of L, fills up R n completely without any overlap, i.e., ∪x∈L {x A common choice of primitive cells is given by i.e., parallelotopes spanned by the primitive vectors of L(A).
The structure element of a crystal is now defined to consist of all points s 1 , . . ., s m that are contained in a particular primitive cell Ξ.Definition 2.3.Let L(A) be a lattice and s ∈ Ξ(A) m , m ∈ N be the structure element.A crystal is defined as the set of tuples The elements of L s (A) are collectively written as x + s = (x + s 1 , x + s 2 , . . ., x + s m ).
Remark 1.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. 1 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
There are infinitely many representations of a crystal L s (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.2, we introduce basic results from integer linear algebra [19], which resolve the relationship of lattice structures. 1 An important tool in this is the notion of a sublattice.
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 The structure of graphene: Illustration of three choices of primitive vectors A, A and A with associated structure elements s, s and s such that L s (A) = L s ( A) = L s ( A).Furthermore, three different primitive cells are illustrated: The corresponding parallelotopes P(A), P( A) as well as a hexagon corresponding to the Voronoi cell of these canonical forms is a crucial ingredient in both the theoretical analysis and the automation of the LFA.
Definition 2.6.A matrix H ∈ R n×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.
Definition 2.7.A matrix D ∈ R n×n is in Smith normal form (SNF) if it is a diagonal matrix and the diagonal entries satisfy d i+1,i+1/d i,i ∈ Z for all i = 1, . . ., n − 1.These entries are called the elementary divisors.
Theorem 2.1.Let A ∈ Q n×n , then (i) there exists a unimodular U ∈ Z n×n such that H = AU is in HNF, (ii) there exist unimodular U,V ∈ Z n×n such that D = V AU is in SNF.
In addition both normal forms H and D are unique with respect to A.
Remark 2. Polynomial algorithms to compute the HNF and SNF can for example be found in [10].
Using these definitions and results one obtains a precise statement about the equality of two lattices.
Theorem 2.2.Let L(A) and L(C) be two lattices then the following statements are equivalent.
(iii) There exists a unimodular matrix U, such that C = AU.
(iv) A and 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.
Definition 2.8.Let L s (A) be a crystal and L(C) ⊂ L(A) be a sublattice.We define the crystal torus T s A,C by . Furthermore, the elements of T s A,C are defined by the equivalence As the primitive vectors of the HNF form a triangular matrix, they allow us to define a canonical ordering of the lattice points on a torus.Corollary 2.1.Let T A,C be arbitrary, i.e., C = AM for some M ∈ Z n×n .Let H ∈ Z n×n be the HNF of M with entries H i j (cf.definition 2.6).Defining the index set I = I 1 × I 2 × . . .× I n by I := {0, 1, . . ., H − 1}, we then obtain In fig. 3 we illustrate corollary 2.1 with an example.Consider Its HNF is given by Thus, a unique list of all representatives of T A,AM is given by In the remainder we drop the bracket notation for reasons of readability.
The H of M yields a lattice basis which allows us to define a canonical lexicographic ordering of the lattice points of a crystal torus T A,AM .

Operators on Crystals
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 with L(Z) ⊂ L(A) being an arbitrary sublattice of L(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 s A,Z , L(Z) ⊂ L(A), precisely describes such an arbitrarily large but finite torus.To shorten notation we use T s A instead of T s A,Z whenever we do not specify Z explicitly.
where T d A corresponds to the crystal of the domain and T c A corresponds to the crystal of the codomain.The function spaces are defined by A value f j (x), x ∈ L(A), corresponds to a value at the position x + s j .The function space is equipped with the scalar product f , g := 1 where f (x), g(x denotes the euclidean scalar product on C |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.[21,Theorem 3.16]), we can connect those operators to the notation of crystal structures.
A ) be a crystal operator.The following statements are equivalent.
1. L is a multiplication operator, i.e., there exist matrices m where the translation operator is defined by For the analysis of such operators the concept of the dual lattice comes in handy as already considered in similar form in [6,26].
A lattice basis of the dual lattice is given by B = A −T .The elements of L(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 section 2. Similar basis functions have already been used in the context of LFA in [3,11].Theorem 3.2.An orthonormal basis for the function space L(T s A,Z ) with a structure element s = (s 1 , . . ., s m ) is given by the wave functions e 1,k , e 2,k , . . ., e m,k defined by e ,k (x Proof.The statement follows by a straightforward, but lengthy calculation, by assuming without loss of generality that M = A −1 Z is given in HNF, making use of corollary 2.1 and the geometric sum formula.
The orthonormal basis of theorem 3.2 can be split into subsets with respect to the wavevector k, i.e., L(T s Theorem 3.3.Let L : L(T s A,Z ) → L(T s A,Z ) be a multiplication operator.Then the subspaces H k of eq.(1) are L-invariant, i.e., L(H k ) ⊆ H k .
Then we obtain by direct calculation Due to their L-invariance the subspaces H k are oftentimes referred to as spaces of harmonics.Thus we can easily represent any A-translationally operator via its symbols, which are formally defined as follows.
) be a multiplication operator with We define the symbol of L according to eq. (2) by In case s = t the spectrum of L can then be extracted from its symbols L k .
Theorem 3.4.Let L : L(T s A,Z ) → L(T s A,Z ) be a multiplication operator with L ∈ C |s|×|s| .
Then spec(L) = ∪ k∈T * A,Z spec(L k ).Proof.Follows immediately due to the orthonormality of the basis e ,k (cf.theorem 3.2) and the L-invariance of the subspaces H k (cf.theorem 3.3).
Remark 3. The main purpose of Z, i.e., the set of primitive vectors that define an arbitrary sublattice L(Z) of L(A), is to simplify the theory developed in section 3 by turning an infinite dimensional setting to an (arbitrarily large) finite one.Though there is another interpretation to it as well.In theorem 3.4 Z explicitly specifies the resolution of the frequency space, i.e., where the spectrum of the multiplication operator is sampled.
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 appendix A 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 2 of multipliers m (y) L −1 = 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 † 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 section 5. 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
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 A,Z with Then, the discretized Laplacian L h : L(T A,Z ) using finite differences is given by The error propagator G of the red-black Gauss-Seidel method can be written with 2 Bounded by the number of lattice points on the (arbitrarily large) torus.respect to the lattice T (0) A,Z ) with where X red and X black correspond to the red and black unknowns of the torus T A,Z as illustrated in fig.4, 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 (0) C,Z = X black and T (0) 4, b)).With respect to this crystal the operators S r and S b can be written as multiplication operators Ŝr , Ŝb : 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., A = C, such that we cannot directly use the computation rules described in appendix A. 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 fully automated as explained in section 5.
As the crystal T s C,Z is yet another representation of T A,Z , the operator L h can be rewritten with respect to this crystal (cf.fig.4) as Lh : Thus, the spectrum of the error propagator of the red-black Gauss-Seidel method applied to the Laplacian Ĝ Ĝk according to eq. ( 3) and the rules in appendix A, followed by the computation of the eigenvalues of the matrices Ĝk .The resulting spectral information of the discretized Laplace operator Lh and the error propagator Ĝ is illustrated in fig.5, where it is sampled on the dual lattice T * C,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
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.
where N B is a diagonal matrix with (N B ) i,i := r • gcd(r, s i ) −1 , where S = V −1 MT −1 = diag (s 1 , . . ., s n ) is the SNF of M (cf.definition 2.7).Consequently, we write L(C) = lcm(L(A), L(B)) and call it the least common multiple of L(A) and L(B).
Proof.Due to lemma 2.1 it is sufficient to find integral matrices N A , N B , such that It can easily be verified that the diagonal matrix N B with the entries (N B ) i,i := r • gcd(r, s i ) −1 guarantees that N A and N B are integral and furthermore minimizes their determinants.Thus, L(C) with C = BT −1 N B is the least common multiple of L(A) and L(B).
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 section 4 for the discretized Laplacian.
defines a structure element of L u (C) such that L u (C) ∼ = L s (A), meant as a one-to-p correspondence.
Proof.Without loss of generality we may assume t j ∈ P(C), j = 1, . . ., p.Then, each element in L s (A) can be written as and there is a unique y, such that Ax = Cy = C y + C(y − y ) with C y ∈ L(C) and C(y − y ) = t j ∈ P(C) ∩ L(A).Thus, we find z as a unique part of the element C y + u = (. . ., C y + t j + s, . ..) = (. . ., z, . ..).
This argument works in the other direction in the same way.
For the sake of clarity of the presentation to follow we opt to write L T s A,C (C) instead of L u (C) as defined in theorem 5.2.Further theorem 5.2 implies a congruence of the function spaces given by the natural isomorphism η : as ( f ) and (η f ) describe the same value distribution on the crystal.In addition this congruence implies that the coarsest possible crystal interpretation is simply the complex coordinate space The scalar product on L(T s A,Z ) then corresponds to the euclidean scalar product on C mn .Using this natural isomorphism of function spaces we can derive similarity transformations for multiplication operators.Theorem 5.3 (Rewriting a multiplication operator with respect to a sublattice).Consider crystals L d (A), L c (A), a sublattice L(C) ⊂ L(A) and a multiplication operator L ∈ C |c|×|d| .
Then, using T A,C = {t 1 , . . ., t p }, the multiplication operator G ∈ C p|c|×p|d| with block matrices (m (y) G Here, the mappings for s ∈ {d, c}, denote the natural isomorphisms between the congruent crystal representations.
Proof.A straightforward calculation for each block-row i yields Due to theorem 5.3 we now know how to rewrite multiple multiplication operators with respect to some common crystal structure with a coarser translational invariance, but it is still not guaranteed that the structure elements are compatible among each other.

A ∼ = T t
A , A ∈ R n are congruent with respect to L(A) if the structure elements are of the same size, i.e., |s| = |t| = m, and there is a permutation π : {1, . . ., m} → {1, . . ., m} as well as shifts y j ∈ L(A), such that 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.

Definition 5.2. Let L : L(T d
A ) → L(T c A ) 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., d i , c j ∈ P(A) = A[0, 1) n for each i, j, • the structure elements d and c are sorted lexicographically. 4e now derive the implications of definition 5.1 for multiplication operators when the structure element is elementwise shifted or permuted.We do so in two steps, theorem 5.4 and theorem 5.5.

Theorem 5.4 (Shifted structure elements). Consider the two multiplication operators
Let further t be a structure element which is obtained from t when shifted element-wise along L(A), i.e., t = (t 1 , . . ., t m ) = ( t1 + y 1 , . . ., tm + y m ) = t + (y 1 , . . ., y m ), where y 1 , . . ., y m ∈ L(A) and m = |t|.Then, the operators L and Ĝ given by fulfill the commutative diagram: Proof.The natural isomorphism between the two corresponding function spaces is given by as f and (T f ) describe the same value distribution on the crystal. 5Again, a straightforward calculation yields Analogously we find Ĝ )g(x + y)] i .Theorem 5.5 (Permuted structure elements).Consider the two multiplication operators L : Let further t be a structure element which is a permuted version of t, i.e., t = ( t1 , . . ., tm ) = (t π(1) , . . ., t π(m) ) = m π t where m = |t|, π : {1, . . ., m} → {1, . . ., m} is a permutation and m π ∈ {0, 1} m×m the corresponding permutation matrix.Then, the operators L and Ĝ given by L ∈ C |t|×|s| with m (y) Proof.Due to the fact that the natural isomorphism p : L(T t A ) → L(T t A ) is a multiplication operator defined by (p f )(x) = m π f (x) for all x ∈ L(A), the statement is true due to the rules of computation in lemma A.1.
Corollary 2.1 and theorems 5.1 to 5.5 allow for the automatic adjustment of crystal representations within the LFA.The corresponding detailed algorithms which make use of these results are given in appendix B.

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 section 2 allow for a canonical description of translationally invariant operators introduced in section 3. Combined with the definition of the dual of a crystal structure in definition 3.2 and a corresponding orthonormal basis in theorem 3.2, the symbols of any single multiplication operator that is translationally invariant with respect to an arbitrary lattice structure, can be expressed (cf. definition 3.3).This combination of choice of basis functions and the explicit connection of operators to their underlying structure enables a completely 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 section 5 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 serve as a tutorial for the use of the algorithms in appendix B.

Multicolored block smoother for the tight-binding Hamiltonian of graphene
In [9] 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 L s (A) where the underlying lattice is triangular, i.e., any three nearby lattice points form an equilateral triangle.We have with the structure element 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. 2. The (nearest neighbor) tight-binding Hamiltonian is a multiplication operator as illustrated in fig.6, a).Overlapping Hexagons We now present an overlapping block smoother for the tightbinding Hamiltonian L. Consider the non-disjoint splitting (or coloring) of the crystal into hexagons as depicted in fig.6, b).This splitting has a translational invariance of C = 2A.Rewriting the graphene crystal T s A with respect to this coarser lattice L(C), we find The splitting is then given by the structure elements 1 , . . ., t 6 ) = (t 2 , t 3 , t 4 , t 5 , t 6 , t 7 ), t (2) = t (1) + a 1 , t (3) = t (1) + a 2 , and t (4) C = and T t (4) C = .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.
A simple way to obtain the corresponding error propagator of the relaxed splitting method which updates the unknowns corresponding to T t ( ) C simultaneously is (in general) as follows.First, we need to find the description L of the underlying operator, the tight-binding Hamiltonian L, with respect to the translational invariance of the splitting 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 (0) L .Consider the structure element t.As can be seen in fig.6, b), the coupling among the unknowns t ( ) which we want to update simultaneously are found in the multipliers: Thus, in order to obtain suitable descriptions for t ( ) , ∈ {2, 3, 4}, we need to consider shifted versions t + τ ( ) .For example τ (1) , τ (2) , τ (3) , τ (4) Finally, the error propagator corresponding to a single color relaxed block Gauss-Seidel update can be written as L m P ) f (x) where m P is the diagonal matrix We summarize the algorithmical steps to obtain the error propagators.

Obtaining error propagators corresponding to splitting methods
Obtain the description of the tight-binding Hamiltonian L with respect to the translational invariance of the splitting via algorithm B.4: Adjust the structure elements, such that the connections among the unknowns updated simultaneously are found within the central multiplier via algorithm B.6: where τ ( ) is defined according to eq. ( 4).Using m P as defined in eq. ( 5) define the operators L( ) m P .
The error propagator of a successive update of the four colors corresponds to the product of the error propagators G ( ) = (I − ω(S ( ) ) † L( ) ).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 algorithm B.1 via ComputeSpectralRadii(g, I, L, S (1) , . . ., S (4) ).
The function g denotes the composition of the error propagators (I, L, S (1) , . . ., S (4) ) 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 section 5.
For ω = 1 2 we obtain the plot in fig.7, a).As the largest spectral radius is greater than one, this method cannot be used as a standalone solver.
Coarse grid correction In [9] a Galerkin coarse grid correction is used with a corresponding coarse grid correction operator that is defined by with L c = RLP and P = R T .Just as in the description of the smoother, the shift invariance of the coarse grid is C = 2A.Defining the coarse structure element by 2s, and with t denoting the structure element of the fine crystal according to the coarse lattice basis, the restriction operator can be described as The multipliers of the restriction operator according to [9] then read The two-grid analysis can now be carried out using algorithm B.1 via ComputeSpectralRadii( f , I, L, S (1) , . . ., S (4) , R).
The function f denotes the composition of the two-grid error propagator (I, L, S (1) , . . ., S (4) , R) A plot of the spectral radii of the two-grid error propagator is given in b) of fig.7 which shows that the two-grid method converges with a convergence rate of ρ max ≈ 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 table 1 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.remark 3).

Two-level analysis for the curl-curl equation
In [3] 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 [7], and a Galerkin coarse grid correction.
The degrees of freedom of the discrete curl-curl equation are associated with the edges of a quadrilateral lattice L( Figure 8: a) Crystal structure of the curl-curl discretization.b) Illustration of the coarse crystal and the restriction operator; dashed lines correspond to restriction weight 1  4 , solid lines to 1  2 .
By multiplying eq. ( 6) with h 2 , the grid size can be incorporated into the material parameter σ h := h 2 σ .Then the discrete operator K = K cc + σ h M can be expressed as a multiplication operator by K : L(T e A ) → L(T e A ) with e = (e h , e v ) = ( where the elements L(A) + e h and L(A) + e v correspond to the (midpoints of the) horizontal and the vertical edge, respectively (cf.fig.8).According to [3] the multipliers are given by: Hybrid smoother The (half) hybrid smoother uses an auxiliary grid.That is, it is a two-grid method itself where smoothing replaces the exact inversion of the auxiliary coarse grid operator.Thus it consists of the following steps: with K c = RKP and P = R T .In here, the coarse crystal is L 2e (2A) and the original crystal L e (A) with respect to the lattice 2A is given by L f (2A) with f = (e, e + a 1 , e + a 2 , e + a 1 + a 2 ).
Then, according to [3] the restriction operator is given as R : L(T f 2A ) → L(T 2s 2A ) with the multipliers: Here, each coarse horizontal and vertical edge is connected to its six nearest edges of the same type (cf.fig.8 b)).The prolongation and coarse grid operator P and K c can be obtained via lemma A.1.
The spectral radii of the two-grid method can then be obtained via where f denotes the composition of the two-grid error propagator Figure 9 b) shows the spectral radii given as sup k ρ(M k ) of the two-grid error propagator as a function of σ h .Here, we actually used the opposite lexicographic orderings in the post-smoothing step to get a result which is similar to [3,Figure 8.1,left].While the results are qualitatively identical, we again attribute the remaining discrepancies to the use of different lexicographic orderings as they have a major impact on the convergence rates.

Conclusion
In this paper we generalized LFA to be 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 translate 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 input is reduced to a bare minimum.Only a representation of the individual operators involved in the analysis of a compounded operator needs to be supplied and each operator can be supplied in its most natural representation.Even though we have limited ourselves in these examples to 2dimensional problems our automated LFA can be applied to operators in higher dimension with ease.In this the only 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.As this cannot be locally resolved in position space, it requires a frequency space mixing analysis and thus does not fit into the presented framework.To be more precise, multiple sequential updates of the same value in one operator lead to linear systems in frequency space that need to be solved to obtain the correct mixing of frequencies (cf.[12,13,20]).While these linear systems could potentially be set up automatically and subsequently be solved by a symbolic calculation, 6 we did not want to include this discussion in this paper.In case this particular calculation would be automated the limitation in the applicability of our approach vanishes.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.section 6.1).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.
We are planning to provide an open source Python implementation of the automated LFA framework in the near future.Annotated Jupyter Notebook excerpts of the examples included in this manuscript can be provided on request.

A Rules of computation
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.(i) Assuming that s = u and t = v, the symbols of L + G are given by L k + G k .
(ii) Assuming that v = s, the symbols of L • G are given by L k • G k .
(iii) The symbols of L * are given by L * k .
(iv) The symbols (L † ) k are given by (L k ) † .

Figure 1 :
Figure 1: Illustration of the crystal representation of three repetitive structures of different origin each consisting of primitive vectors a 1 , a 2 , a shaded primitive cell P(A) and a structure element (s 1 , s 2 ).a) Staggered discretization of the curl-curl system of PDEs with unknowns on the edges of a rectangular lattice.b) Hexagonal lattice of the carbon allotrope graphene.c) Checkerboard coloring of a rectangular lattice as it is encountered in the red-black Gauss-Seidel smoother.

Figure 4 :
Figure 4: The discretized Laplace operator with respect to two different crystal representations.In a) the operator is illustrated with respect to the primitive vectors a 1 , a 2 and in b) with respect to c 1 , c 2 .

Theorem 5 . 1 .
Given two n-dimensional Lattices L(A), L(B).If there exists an r ∈ Z, such that M = rA −1 B ∈ Z n×n , then there is a smallest lattice L(C) with L(C) ⊂ L(A) and L(C) ⊂ L(B).A lattice basis of L(C) is given by

Figure 5 :
Figure 5: The spectra of a) the discretized Laplace operator and b) the red-black Gauss-Seidel smoother with respect to the red-black crystal structure.In here, b 1 and b 2 denote the columns of C −T .

Figure 6 :
Figure 6: a) Schematic stencil of the tight binding Hamiltonian L [t 0 ,t 1 ] of graphene.b) Illustration of the domain decomposition into hexagons.Each unknown belongs to 3 different colors.

4 .
. Smooth on K N x N = r N with zero initial guess: x N ← S † N r N , Prolongate to the primal crystal x ← x + P N x N .In fig.9 a) a contourplot of ρ(G k ) with σ h = 0.01, with respect to k ∈ P(A −T ) is given.It corresponds to the result given in [3, Figure6.2, right].These plots differ only by a reflection along the main diagonal which we suspect is a result from the use of different lexicographic orderings.Coarse grid correction In[3] a Galerkin coarse grid correction is used corresponding to the error propagator E = (I − P(K c ) † RK)

Figure 9 :
Figure 9: a) Contourplot of ρ(G k ) for the hybrid smoother with respect to k ∈ A T = A, σ h = 0.01.b) Contourplot of the convergence estimates of the two-grid method sup k ρ(M k ) with respect to σ h .The plots correspond to the results [3, Figure 6.2, right] and [3, Figure 8.1, left], respectively.

Lemma A. 1 . 1 .G
Let two multiplication operators be given byL : L(T s A,Z ) → L(T t A,Z ), (L f )(x) = ∑ y∈T A,Z m (y) L f (x + y), m(y)L ∈ C |t|×|s| , G : L(T u A,Z ) → L(T v A,Z ), (G f )(x) = ∑ y∈T A,Z m (y) L f (x + y), m(y)G ∈ C |v|×|u| .Then the following operators are multiplication operators as well:(i) If s = u and t = v, then L + G : L(T s A,Z ) → L(T t A,Z ) with m If v = s, then LG : L(T u A,Z ) → L(T t A,Z) with m The adjoint is given by L * : L(T t A,Z ) → L(T s A,Z ) with m Let two multiplication operators be given byL : L(T s A,Z ) → L(T t A,Z ), (L f )(x) = ∑ y∈T A,Z m (y) L f (x + y), m(y)L ∈ C |t|×|s| , G : L(T u A,Z ) → L(T v A,Z ), (G f )(x) = ∑ ∈ C |v|×|u|with corresponding symbols L k and G k .Then we have the following statements.

Table 1 :
ρ analytic | Convergence history of the two-grid method applied to the tight-binding Hamiltonian of graphene with (41 × 41) • 23unknowns/atoms and periodic boundary conditions.The reported asymptotic convergence rate ρ i coincides with high precision to the convergence estimate ρ analytic = sup k∈T *