Block Generalized Locally Toeplitz Sequences : From the Theory to the Applications

The theory of generalized locally Toeplitz (GLT) sequences is a powerful apparatus for computing the asymptotic spectral distribution of matrices An arising from virtually any kind of numerical discretization of differential equations (DEs). Indeed, when the mesh fineness parameter n tends to infinity, these matrices An give rise to a sequence {An}n, which often turns out to be a GLT sequence or one of its “relatives”, i.e., a block GLT sequence or a reduced GLT sequence. In particular, block GLT sequences are encountered in the discretization of systems of DEs as well as in the higher-order finite element or discontinuous Galerkin approximation of scalar DEs. Despite the applicative interest, a solid theory of block GLT sequences has been developed only recently, in 2018. The purpose of the present paper is to illustrate the potential of this theory by presenting a few noteworthy examples of applications in the context of DE discretizations.


Introduction
The theory of generalized locally Toeplitz (GLT) sequences stems from Tilli's work on locally Toeplitz (LT) sequences [39] and from the spectral theory of Toeplitz matrices [2,[10][11][12][13]30,32,38,40,41,43].It was then carried forward in [26,27,36,37], and it has been recently extended by Barbarino in [3].This theory is a powerful apparatus for computing the asymptotic spectral distribution of matrices arising from the numerical discretization of continuous problems, such as integral equations (IEs) and, especially, differential equations (DEs).The experience reveals that virtually any kind of numerical methods for the discretization of DEs gives rise to structured matrices A n whose asymptotic spectral distribution, as the mesh fineness parameter n tends to infinity, can usually be computed through the theory of GLT sequences.We refer the reader to [26,Section 10.5], [27,Section 8.5], and [9,36,37] for applications of the theory of GLT sequences in the context of finite difference (FD) discretizations of DEs; to [26,Section 10.6] and [5,9,37] for the finite element (FE) case; to [7] for the finite volume (FV) case; to [26,Section 10.7], [27,Section 8.5], and [17,22,24,25,33] for the case of isogeometric analysis (IgA) discretizations, both in the collocation and Galerkin frameworks; and to [18] for a further recent application to fractional DEs.We also refer the reader to [26,Section 10.4] and [1,34] for a look at the GLT approach for sequences of matrices coming from IE discretizations.
In the very recent work [29], starting from the original intuition by the third author [37,Section 3.3], the theory of block GLT sequences has been developed in a systematic way as an extension of the theory of GLT sequences.Such an extension is of the utmost importance in practical applications.In particular, it provides the necessary tools for computing the spectral distribution of block structured matrices arising from the discretization of systems of DEs [37,Section 3.3] and from the higher-order finite element or discontinuous Galerkin approximation of scalar/vectorial DEs [6,20,28].The purpose of this paper is to illustrate the potential of the theory of block GLT sequences [29] and of its multivariate version -which combines the results of [29] with the "multivariate technicalities" from [27] -by presenting a few noteworthy examples of applications.Actually, the present paper can be seen as a necessary completion of the purely theoretical work [29].
The paper is organized as follows.In Section 2 we report a summary of the theory of block GLT sequences.In Section 3 we focus on the FD discretization of a model system of univariate DEs; through the theory of block GLT sequences, we compute the spectral distribution of the related discretization matrices.In Section 4 we focus on the higher-order FE approximation of the univariate diffusion equation; again, we compute the spectral distribution of the associated discretization matrices through the theory of block GLT sequences.In Section 5 we summarize the multivariate version of the theory of block GLT sequences, also known as the theory of multilevel block GLT sequences.In Section 6 we describe the general GLT approach for computing the spectral distribution of matrices arising from the discretization of systems of partial differential equations (PDEs).In Section 7 we focus on the B-spline IgA approximation of a bivariate variational problem for the curl-curl operator, which is of interest in magnetostatics; through the theory of multilevel block GLT sequences, we compute the spectral distribution of the related discretization matrices.In Section 8 we focus on the linear elasticity equations discretized through the modified FE Taylor-Hood Q1isoQ1 stable pair; again, we compute the spectral distribution of the related discretization matrices through the theory of multilevel block GLT sequences.Final considerations are collected in Section 9.

The Theory of Block GLT Sequences
In this section we summarize the theory of block GLT sequences, which was originally introduced in [37, Section 3.3] and has been recently revised and systematically developed in [29].
Sequences of Matrices and Block Matrix-Sequences.Throughout this paper, a sequence of matrices is any sequence of the form {A n } n , where A n is a square matrix of size d n and d n → ∞ as n → ∞.Let s ≥ 1 be a fixed positive integer independent of n; an s-block matrix-sequence (or simply a matrix-sequence if s can be inferred from the context or we do not need/want to specify it) is a special sequence of matrices {A n } n in which the size of A n is d n = sn.

Singular Value and Eigenvalue Distribution of a Sequence of Matrices.
Let µ k be the Lebesgue measure in R k .Throughout this paper, all the terminology from measure theory (such as "measurable set", "measurable function", "a.e.", etc.) is referred to the Lebesgue measure.A matrix-valued function f : D ⊆ R k → C r×r is said to be measurable (resp., continuous, Riemann-integrable, in L p (D), etc.) if its components f αβ : D → C, α, β = 1, . . ., r, are measurable (resp., continuous, Riemann-integrable, in L p (D), etc.).We denote by C c (R) (resp., C c (C)) the space of continuous complex-valued functions with bounded support defined on R (resp., C).If A ∈ C m×m , the singular values and the eigenvalues of A are denoted by σ 1 (A), . . ., σ m (A) and λ 1 (A), . . ., λ m (A), respectively.Definition 1.Let {A n } n be a sequence of matrices, with A n of size d n , and let f : D ⊂ R k → C r×r be a measurable function defined on a set D with 0 < µ k (D) < ∞.
• We say that {A n } n has a (asymptotic) singular value distribution described by f , and we write {A n } n ∼ σ f , if • We say that {A n } n has a (asymptotic) spectral (or eigenvalue) distribution described by f , and we write If {A n } n has both a singular value and an eigenvalue distribution described by f , we write {A n } n ∼ σ,λ f .
We note that Definition 1 is well-posed because the functions taking values in C r×r for some r ≥ 1.We refer the reader to [28,Remark 1] or [21,Section 4] for the informal meaning behind the spectral distribution (2); a completely analogous meaning can also be given for the singular value distribution (1).
The next two theorems are useful tools for computing the spectral distribution of sequences formed by Hermitian or perturbed Hermitian matrices.For the related proofs, we refer the reader to [31,Theorem 4.3] and [4,Theorem 1.1].In what follows, the conjugate transpose of the matrix A is denoted by A * .If A ∈ C m×m and 1 ≤ p ≤ ∞, we denote by A p the Schatten p-norm of A, i.e., the p-norm of the vector (σ 1 (A), . . ., σ m (A)).The Schatten ∞-norm A ∞ is the largest singular value of A and coincides with the spectral norm A .The Schatten 1-norm A 1 is the sum of the singular values of A and is often referred to as the trace-norm of A. The Schatten 2-norm A 2 coincides with the Frobenius norm of A. For more on Schatten p-norms, see [8].• the matrices X n are Hermitian and ), its Fourier coefficients are denoted by where the integrals are computed componentwise.The nth block Toeplitz matrix generated by f is defined as It is not difficult to see that all the matrices T n ( f ) are Hermitian when f is Hermitian a.e.
Block Diagonal Sampling Matrices.For n ∈ N and a : [0, 1] → C s×s , we define the block diagonal sampling matrix D n (a) as the diagonal matrix Zero-Distributed Sequences.A sequence of matrices {Z n } n such that {Z n } n ∼ σ 0 is referred to as a zero-distributed sequence.Note that, for any r ≥ 1, (throughout this paper, O m and I m denote the m × m zero matrix and the m × m identity matrix, respectively).Proposition 1 provides an important characterization of zero-distributed sequences together with a useful sufficient condition for detecting such sequences.Throughout this paper we use the natural convention 1/∞ = 0.
Proposition 1.Let {Z n } n be a sequence of matrices, with Z n of size d n .
Approximating Classes of Sequences.The notion of approximating classes of sequences (a.c.s.) is the fundamental concept on which the theory of block GLT sequences is based.
Definition 2. Let {A n } n be a sequence of matrices, with A n of size d n , and let {{B n,m } n } m be a sequence of sequences of matrices, with B n,m of size d n .We say that {{B n,m } n } m is an approximating class of sequences (a.c.s.) for {A n } n if the following condition is met: for every m there exists n m such that, for n ≥ n m , where n m , c(m), ω(m) depend only on m and lim Roughly speaking, {{B n,m } n } m is an a.c.s. for {A n } n if, for large m, the sequence {B n,m } n approximates {A n } n in the sense that A n is eventually equal to B n,m plus a small-rank matrix (with respect to the matrix size d n ) plus a small-norm matrix.It turns out that, for each fixed sequence of positive integers d n such that d n → ∞, the notion of a.c.s. is a notion of convergence in the space If X ∈ C m 1 ×m 2 and Y ∈ C 1 × 2 are any two matrices, the tensor (Kronecker) product of X and Y is the m 1 1 × m 2 2 matrix defined as follows: We recall that the tensor product operation ⊗ is associative and bilinear.Moreover, Finally, if X 1 , X 2 can be multiplied and Y 1 , Y 2 can be multiplied, then Lemma 1.For i, j = 1, . . ., s, let {A n,ij } n be a sequence of matrices and suppose that {B Proof.Let E ij be the s × s matrix having 1 in position (i, j) and 0 elsewhere.Note that Since {B −→ {A n,ij } n , it is clear from (3), (4) and the definition of a.c.s. that Now, if {B n } n (this an obvious consequence of the definition of a.c.s.).Thus, the thesis follows from (7) and (8).
Block GLT Sequences.Let s ≥ 1 be a fixed positive integer.An s-block GLT sequence (or simply a GLT sequence if s can be inferred from the context or we do not need/want to specify it) is a special s-block matrix-sequence {A n } n equipped with a measurable function κ : [0, 1] × [−π, π] → C s×s , the so-called symbol.We use the notation The main properties of s-block GLT sequences proved in [29] are listed below.If A is a matrix, we denote by A † the Moore-Penrose pseudoinverse of A (recall that A † = A −1 whenever A is invertible).If f m , f : D ⊆ R k → C r×r are measurable matrix-valued functions, we say that f m converges to f in measure (resp., a.e., in L p (D), etc.) if ( f m ) αβ converges to f αβ in measure (resp., a.e., in L p (D), etc.) for all α, β = 1, . . ., r. −→ {A n } n and κ m → κ in measure.

FD Discretization of a System of DEs
Consider the following system of DEs: x ∈ (0, 1), x ∈ (0, 1), In this section we consider the classical central FD discretization of (9).Through the theory of block GLT sequences we show that the corresponding sequence of (normalized) FD discretization matrices enjoys a spectral distribution described by a 2 × 2 matrix-valued function, where 2 is the number of equations that compose the system (9).In what follows, we use the following notation: tridiag j=1,...,n The parameters α j , β j , γ j may be either scalars or s × s blocks for some s > 1, in which case the previous matrix is a block tridiagonal matrix.
• The unknowns are sorted as follows: • The equations are sorted as follows, in accordance with the ordering (15) for the unknowns: Suppose we decide to change the ordering for both the unknowns and the equations.More precisely, suppose we opt for the following orderings.
• The unknowns are sorted as follows: • The equations are sorted as follows, in accordance with the ordering (17) for the unknowns: The matrix C n associated with the linear system (10) assuming the new orderings ( 17) and ( 18) is the 2 × 2 block tridiagonal matrix given by The matrix C n is similar to B n .Indeed, by permuting both rows and columns of B n according to the permutation 1, n + 1, 2, n + 2, . . ., n, 2n we obtain C n .More precisely, let e 1 , . . ., e n and ẽ1 , . . ., ẽ2n be the vectors of the canonical basis of R n and R 2n , respectively, and let Π n be the permutation matrix associated with the permutation 1, n + 1, 2, n + 2, . . ., n, 2n, that is, . . Then
Theorem 3. Suppose that a 11 , a 12 , a 21 , and If moreover a 21 = −a 12 then we also have Proof.From (19) we have where E pq is the 2 × 2 matrix having 1 in position (p, q) and 0 elsewhere.It is clear that, for every As a consequence, the decomposition (28), GLT 2 and GLT 3 imply (25), which in turn implies (26) by GLT 1.It only remains to prove (27) in the case where a 21 = −a 12 .In this case, we have Consider the symmetric approximation of C n given by see, e.g., [26,Section 2.4.1].Therefore: Thus, (27) follows from Theorem 2.

Higher-Order FE Discretization of the Diffusion Equation
Consider the diffusion equation In this section we consider the higher-order FE discretization of (30).Through the theory of block GLT sequences we show that the corresponding sequence of (normalized) FE discretization matrices enjoys a spectral distribution described by a (p − k) × (p − k) matrix-valued function, where p and k represent, respectively, the degree and the smoothness of the piecewise polynomial functions involved in the FE approximation.Note that this result represents a remarkable argument in support of [28, Conjecture 2].

FE Discretization
The weak form of (30) reads as follows: In the FE method we fix a set of basis functions {ϕ 1 , . . ., ϕ N } ⊂ H 1 0 ([0, 1]) and we look for an approximation of the exact solution in the space W = span(ϕ 1 , . . ., ϕ N ) by solving the following discrete problem: find Since {ϕ 1 , . . ., ϕ N } is a basis of W , we can write u W = ∑ N j=1 u j ϕ j for a unique vector u = (u 1 , . . ., u N ) T .By linearity, the computation of u W (i.e., of u) reduces to solving the linear system T and A is the stiffness matrix,

p-Degree C k B-spline Basis Functions
Following the higher-order FE approach, the basis functions ϕ 1 , . . ., ϕ N will be chosen as piecewise polynomials of degree p ≥ 1.More precisely, for p, n ≥ 1 and 0 k] : R → R be the B-splines of degree p and smoothness C k defined on the knot sequence We collect here a few properties of B 1,[p,k] , . . ., B n(p−k)+k+1,[p,k] that we shall use in this paper.For the formal definition of B-splines, as well as for the proof of the properties listed below, see [16,35].
• The support of the ith B-spline is given by • Except for the first and the last one, all the other B-splines vanish on the boundary of [0, 1], i.e., • {B 1,[p,k] , . . . ,B n(p−k)+k+1,[p,k] } is a basis for the space of piecewise polynomial functions of degree p and smoothness C k , that is, n ] ∈ P p for all i = 0, . . ., n − 1 , where P p is the space of polynomials of degree ≤ p.Moreover, {B 2,[p,k] , . . ., B n(p−k)+k,[p,k] } is a basis for the space • The B-splines form a non-negative partition of unity over [0, 1]: • The derivatives of the B-splines satisfy where c p is a constant depending only on p.Note that the derivatives B i,[p,k] may not be defined at some of the grid points 0, 1 n , 2 n , . . ., n−1 n , 1 in the case of C 0 smoothness (k = 0).In (37) it is assumed that the undefined values are excluded from the summation.In formulas, setting for the B-splines We point out that the supports of the reference B-splines The basis functions ϕ 1 , . . ., ϕ N are defined as follows: In particular, with the notations of Section 4.

GLT Analysis of the Higher-Order FE Discretization Matrices
The stiffness matrix (31) resulting from the choice of the basis functions as in (40) will be denoted by A n,[p,k] (a), The main result of this section (Theorem 4) gives the spectral distribution of the normalized sequence {n −1 A n,[p,k] (a)} n .The proof of Theorem 4 is entirely based on the theory of block GLT sequences and it is therefore referred to as "GLT analysis".It also requires the following lemma, which provides an approximate construction of the matrix A n,[p,k] (1) corresponding to the constant-coefficient case where a(x) = 1 identically.In view of what follows, define the (p and the (p Due to the compact support of the reference functions β Proof.By ( 33) and ( 39), for all r, R = 1, . . ., n − ν and q, Q = 1, . . ., p − k we have which completes the proof.
Proof.The proof consists of four steps.Throughout this proof, we will use the following notations.
• For every square matrix A of size n(p − k) + k − 1 we denote by A the principal submatrix of A corresponding to the row and column indices i, j Step 1.Consider the linear operator .
The next three steps are devoted to show that Once this is done, the theorem is proved.Indeed, from (44) we immediately obtain the relation Step 2. We first prove (44) in the constant-coefficient case where g(x) = 1 identically.In this case, by Lemma 2 we have n Step 3. Now we prove (44) in the case where g ∈ C([0, 1]).Let By (32), ( 33) and ( 37), for all r, R = 1, . . ., n − ν and q, Q = 1, . . ., p − k we have where ω g (•) is the modulus of continuity of g and the last inequality is justified by the fact that the distance of the point r/(n − ν) from the interval [(r − Step 4. Finally, we prove (44) in the general case where g ∈ L 1 ([0, 1]).By the density of Moreover, We show that {n Once this is done, the thesis (44) follows immediately from GLT 4. To prove (47), we recall that see, e.g., [26,Section 2.4.3].By (37) we obtain Thus, the a.c.s.convergence (47) follows from Proposition 2.
Remark 1.By following step by step the proof of Theorem 4, we can give an alternative (much simpler) proof of [6, Theorem A.6] based on the theory of block GLT sequences.

The Theory of Multilevel Block GLT Sequences
As illustrated in Sections 3 and 4, the theory of block GLT sequences allows the computation of the singular value and eigenvalue distribution of block structured matrices arising from the discretization of univariate DEs.In order to cope with multivariate DEs, i.e., PDEs, we need the multivariate version of the theory of block GLT sequences, also known as the theory of multilevel block GLT sequences.The present section is devoted to a careful presentation of this theory, which is obtained by combining the results of [29] with the necessary technicalities for tackling multidimensional problems [27].
Multi-Index Notation.The multi-index notation is an essential tool for dealing with sequences of matrices arising from the discretization of PDEs.A multi-index i ∈ Z d , also called a d-index, is simply a (row) vector in Z d ; its components are denoted by i 1 , . . ., i d .
• For any d-index m, we set N(m) = ∏ d j=1 m j and we write m → ∞ to indicate that min(m) → ∞.
We assume for this set the standard lexicographic ordering: . . .[ (j 1 , . . ., For instance, in the case d = 2 the ordering is the following: • When a d-index j varies over a multi-index range h, . . ., k (this is sometimes written as j = h, . . ., k), it is understood that j varies from h to k following the specific ordering (49).For instance, if m ∈ N d and if we write x = [x i ] m i=1 , then x is a vector of size N(m) whose components x i , i = 1, . . ., m, are ordered in accordance with (49): the first component is x 1 = x (1,...,1,1) , the second component is x (1,...,1,2) , and so on until the last component, which is , then X is a N(m) × N(m) matrix whose components are indexed by two d-indices i, j, both varying from 1 to m according to the lexicographic ordering (49).
• Given h, k ∈ Z d with h ≤ k, the notation ∑ k j=h indicates the summation over all j in h, . . ., k. • Operations involving d-indices that have no meaning in the vector space Z d must be interpreted in the componentwise sense.For instance, ij = (i 1 j 1 , . . ., i d j d ), i/j = (i 1 /j 1 , . . ., i d /j d ), etc.
Multilevel Block Matrix-Sequences.Given d, s ≥ 1, a d-level s-block matrix-sequence (or simply a matrix-sequence if d and s can be inferred from the context or we do not need/want to specify them) is a sequence of matrices of the form {A n } n , where: • n varies in some infinite subset of N; , its Fourier coefficients are denoted by where k • θ = k 1 θ 1 + . . .+ k d θ d and the integrals are computed componentwise.For n ∈ N d , the nth multilevel block Toeplitz matrix generated by f is defined as It is not difficult to see that the map f → T n ( f ) is linear.Moreover, it can be shown that where the transpose conjugate function f * is defined by f * (θ) = ( f (θ)) * ; in particular, all the matrices T n ( f ) are Hermitian whenever f is Hermitian a.e.We also recall that if n ∈ N d and f 1 , f 2 , . . ., where The main properties of d-level s-block GLT sequences are listed below.
We have:

Discretizations of Systems of PDEs: The General GLT Approach
In this section we outline the main ideas of a multidimensional block GLT analysis for general discretizations of PDE systems.What we are going to present here is then a generalization of what we have seen in Section 3. We begin by proving a series of auxiliary results.In what follows, given n ∈ N d and s ≥ 1, we denote by Π n,s the permutation matrix given by where e i , i = 1, . . ., n, are the vectors of the canonical basis of R N(n) , which, for convenience, are indexed by a d-index i = 1, . . ., n instead of a linear index i = 1, . . ., N(n).Note that Π n,2 coincides with the matrix Π n in (20).
is similar via the permutation (52) to the multilevel block Toeplitz matrix Proof.Let E ij be the s × s matrix having 1 in position (i, j) and 0 elsewhere.Since By ( 5) and ( 6), Proof.With obvious adaptations, it is the same as the proof of Lemma 3.
We recall that a d-variate trigonometric polynomial is a finite linear combination of the d-variate Fourier frequencies e ik•θ , k ∈ Z d .Theorem 5.For i, j = 1, . . ., s, let {A n,ij } n be a d-level 1-block GLT sequence with symbol . Then, the matrix-sequence {Π n,s A n Π T n,s } n is a d-level s-block GLT sequence with symbol κ.
Proof.The proof consists of the following two steps.
Step 1.We first prove the theorem under the additional assumption that A n,ij is of the form where .
By setting L = max i,j=1,...,s L ij and by adding zero matrices of the form D n (0)T n (0) in the summation (53) whenever L ij < L, we can assume, without loss of generality, that with L independent of i, j.Let E ij be the s × s matrix having 1 in position (i, j) and 0 elsewhere.Then, By Lemmas 3 and 4, It follows that Thus, by GLT 2 and GLT 3, {Π n,s A n Π T n,s } n is a d-level s-block GLT sequence with symbol Step 2. We now prove the theorem in its full generality.Since {A n,ij } n ∼ GLT κ ij , by [27,Corollary 7.2] there exist functions a ,ij , = 1, . . ., ,ij is a d-variate trigonometric polynomial, • κ . We have: We conclude that {Π n,s A n Π T n,s } n ∼ GLT κ by GLT 4.
Now, suppose we have a system of linear PDEs of the form where x ∈ (0, Moreover, it turns out that, after a suitable normalization that we ignore in this discussion, 1 A n has the following block structure: where A n,ij is the (normalized) matrix coming from the discretization of the differential operator L ij .
In the simplest case where the operators L ij have constant coefficients and we use equispaced grids in each direction, the matrix A n,ij takes the form where f ij is a d-variate trigonometric polynomial, while the perturbation Z n,ij is usually a low-rank correction due to boundary conditions and, in any case, we have {Z n,ij } n ∼ σ 0. Hence, by GLT 2 and GLT 3, and it follows from Theorem 5 that In the case where the operators L ij have variable coefficients, the matrix A n,ij usually takes the form where by GLT 2 and GLT 3, and it follows from Theorem 5 that 1 The normalization we are talking about is the analog of the normalization that we have seen in Section 3, which allowed us to pass from the matrix A n in (12) to the matrix B n in (14). 2 For example, in Section 3, while proving (21), we have seen that K n (a 11 ), which plays there the same role as the matrix A n,11 here, is equal to D n (a 11 )T n (2 − 2 cos θ) + Z n for some zero-distributed sequence {Z n } n .

Compatible B-Spline IgA Discretization
Let n = (n 1 , n 2 ) ∈ N 2 , let p ≥ 2, and define the space Following a compatible B-spline approach [15], we look for an approximation of the solution in the space V n,[p] (curl, Ω) by solving the following discrete problem: find After choosing a suitable ordering on the basis functions of V n,[p] (curl, Ω) displayed in (61), by linearity the computation of u V reduces to solving a linear system whose coefficient matrix is given by where Note that M n,[p−1] is a square matrix of size n + p − 1, K n,[p] is a square matrix of size n + p, while H n,[p] is a rectangular matrix of size (n + p − 1) × (n + p).

GLT Analysis of the B-Spline IgA Discretization Matrices
In the main result of this section (Theorem 6), assuming that n = nν for a fixed vector ν, we show that the spectral distribution of the sequence {A n,[p] } n is described by a 2 × 2 matrix-valued function whose determinant is zero everywhere (Remark 2).In order to prove Theorem 6, some preliminary work is necessary.We first note that, in view of the application of Theorem 5, the matrix A n,[p] has an unpleasant feature: the anti-diagonal blocks A n,[p], 12 and A n,[p where M n,[p−1] and H n,[p] are square matrices of size n + p given by In view of the application of Theorem 1, we note that where and we note the three series are actually finite sums because of (57).

Linear Elasticity Equations Discretized by the Modified Taylor-Hood FE Q1isoQ1 Stable Pair
In this section we focus on a specific formulation of the linear elasticity equations, which is related to the glacial isostatic adjustment (GIA) model; see, e.g., [42].The GIA model is used to describe the response of the Earth to the redistribution of mass due to alternating glaciation and deglaciation periods.The considered elasticity problem is the following bidimensional system of PDEs, which is defined over a rectangular domain Ω ⊂ R 2 and is accompanied by appropriate boundary conditions.In (77), u = [u 1 (x 1 , x 2 ), u 2 (x 1 , x 2 )] T is a vector containing the displacements u 1 and u 2 in the x 1 and x 2 directions, respectively; p is the pressure; b = [b 1 , b 2 ] T is a constant advection vector; and ρ = α 2 /λ, where α and λ are the so-called Lamé coefficients, which are assumed to be constant.
Note that problem (77) is of the form (54). Hence, we expect that a standard discretization of (77) yields a matrix having a block structure of the form (55). Indeed, if we discretize (77) by the modified Taylor-Hood Q1isoQ1 FE stable pair [14, p. 167] on a uniform mesh with stepsize 1 n in each direction x i , the computation of the numerical solution reduces to solving a linear system whose coefficient matrix is a two-by-two block matrix given by where n = (n, n).Here, the pivot block K n is a two-by-two block matrix, M n is the so-called mass matrix, and the blocks } n ∼ σ f nonsym (θ 1 , θ 2 ).
As a consequence, thanks to GLT 3, If we include the boundary effects, the Q1isoQ1 FE scheme for discretizing the original problem yields a K n which, up to a permutation, is given by T n ( f K ) + E n , where E n is a low-rank perturbation.Thus, by Theorem 5, the symbol of {K n } n remains the same, regardless of the boundary conditions.As for the matrix K n , also in this case the symbol does not change when we include the boundary conditions.
Symbol of the anti-diagonal part of A n .Due to the rectangular nature of the blocks H n , H T n , we cannot compute their symbols directly.Our strategy is therefore to embed them into larger matrices and to derive their symbols by applying Theorem 1 on the anti-diagonal part of A n , which is denoted by J n in what follows.We begin by rewriting H n,1 , H n,2 as the downsampling of larger square matrices H n,1 , H n,2 , namely where H n,k are five-diagonal block matrices such that each block is itself five-diagonal, while C is called cutting matrix and is a special projection matrix that deletes columns block-wise and within the blocks of H n,k .More in detail, n H n,1 , n H n,2 are 2-level Toeplitz matrices associated with the following bivariate complex-valued functions: Therefore, if we denote by J n the anti-diagonal two-by-two block matrix with J n,21 = H n , J n,12 = H T n , and H n = [ H n,1 , H n,2 ], and we use relation (78) combined with Theorem 5 and GLT 1, we obtain As a consequence, by observing that and by applying Theorem 1, we arrive at {nJ n } n ∼ σ,λ f J (θ 1 , θ 2 ).

Conclusions
We have illustrated through specific examples the application potential of the theory of block GLT sequences and of its multivariate version, thus bringing to completion the purely theoretical work [29].It should be said, however, that the theory of GLT sequences is still incomplete.In particular, besides filling in the details of the theory of multilevel block GLT sequences, 3 it will be necessary to develop the theory of the so-called reduced GLT sequences, as explained in [26,Chapter 11].
as the principal submatrix corresponding to the row and column indices i, j = k + 1, . . ., k + (n − ν)(p − k) and zeros elsewhere.Note that P T n,[p,k] P n,[p,k] = I (n−ν)(p−k) and P T n,[p,k] AP n,[p,k] = A for every square matrix A of size n(p − k) + k − 1.

Lemma 4 .
Let n ∈ N d , let a ij : [0, 1] d → C for i, j = 1, . . ., s, and set a = [a ij ] s i,j=1 .The block matrix D n = [D n (a ij )] s i,j=1 is similar via the permutation (52) to the multilevel block diagonal sampling matrix D n (a), that is, Π n,s D n Π T n,s = D n (a).
e. (and hence also in measure);• {Π n,s A (m) n Π T n,s } n a.c.s.−→ {Π n,s A n Π T n,s } n because {A (m) n } n a.c.s.−→ {A n } n by Lemma 1.

Theorem 1 .
Let{X n } n be a sequence of matrices, with X n Hermitian of size d n , and let {P n } n be a sequence such that P n ∈ C d n ×δ n , P * n P n = I δ n , δ n ≤ d n and δ n /d n → 1 as n → ∞.Then, {X n } n ∼ σ,λ κ if and only if {P * n X n P n } n ∼ σ,λ κ.Theorem 2. Let {X n } n and {Y n } n be sequences of matrices, with X n and Y n of size d n .Assume that: [26, precisely, there exists a pseudometricd a.c.s. in E such that {{B n,m } n } m is an a.c.s. for {A n } n if and only if d a.c.s.({B n,m } n , {A n } n ) → 0 as m → ∞.We will therefore use the convergence notation {B n,m } n → {A n } n to indicate that {{B n,m } n } m is an a.c.s. for {A n } n .A useful criterion to identify an a.c.s. is provided in the next proposition[26, Corollary 5.3]. − that κ is invertible a.e.GLT 4. {A n } n ∼ GLT κ if and only if there exist s-block GLT sequences {B n,m } n ∼ GLT κ m such that {B n,m } n Figure 1.B-splines B 1,[p,k] , . . ., B n(p−k)+k+1,[p,k] for p = 3 and k = 1, with n = 10.
(29)ollows that each entry of Z n,[p,k] (g) is bounded in modulus by C p ω g (1/n), where C p is a constant depending only on p.Moreover, by(33), the matrix Z n,[p,k] (g) is banded with bandwidth bounded by a constant w p depending only on p.Thus, by(29), Z n,[p,k] (g) ≤ w p C p ω g (1/n) → 0 as n → ∞, and so {Z n,[p,k] (g)} n is zero-distributed by Proposition 1.Since [27, e.g.,[27, Lemma 5.3].For n ∈ N d and a : [0, 1] d → C s×s , we define the multilevel block diagonal sampling matrix D n (a) as the block diagonal matrix Multilevel Block GLT Sequences.Let d, s ≥ 1 be fixed positive integers.A d-level s-block GLT sequence (or simply a GLT sequence if d and s can be inferred from the context or we do not need/want to specify them) is a special d-level s-block matrix-sequence {A n } n equipped with a measurable function κ : 1) d .The matrices A n resulting from any standard discretization of (54) are parameterized by a d-index n = (n 1 , . . ., n d ), where n i is related to the discretization step h i in the ith direction, and n i → ∞ if and only if h i → 0 (usually, h i ≈ 1/n i ).By choosing each n i as a function of a unique discretization parameter n ∈ N, as it normally happens in practice where the most natural choice is n i = n for all i = 1, . . ., d, we see that n = n(n) and, consequently, {A n } n is a (d-level) matrix-sequence.