Domain decomposition methods are well-studied approaches for the numerical solution of partial differential equations both experimentally and theoretically[5, 13, 15, 36], due to their efficiency and robustness for many large-scale problems, and the need for parallel algorithms. Among the many families of domain decomposition algorithms are Neumann-Neumann , FETI-DP , Schwarz [15, 36], and Optimized Schwarz [13, 20]. Balancing domain decomposition by constraints (BDDC) is one family of non-overlapping domain decomposition methods. While BDDC was first introduced by Dohrmann in , several variants have recently been proposed. BDDC-like methods have been successfully applied to many PDEs, including elliptic problems , the incompressible Stokes equations [24, 26], H(curl) problems , flow in porous media , and the incompressible elasticity problem [11, 32]. Theoretical analysis of BDDC has primarily been based on finite-element approximation theory [8, 14, 29, 30]. It has been shown that the condition number of the preconditioned BDDC operator can be bounded by a function of (where is the meshsize, and is the subdomain size), independent of the number of subdomains [29, 30]. A nonoverlapping domain decomposition method for discontinuous Galerkin based on the BDDC algorithm is presented in 
, and the condition number of the preconditioned system is shown to be bounded by similar estimates as those for conforming finite element methods. BDDC methods in three- or multilevel forms have also been developed[31, 39, 40], and good implementations are available, for example, [1, 34, 43].
Since BDDC algorithms are widely used to solve many problems with high efficiency and parallelism, better understanding of how this methodology works is useful in the design of new algorithms. Local Fourier analysis (LFA), first introduced by Brandt  and well-studied for multigrid methods [9, 35, 37, 41, 42], is an analysis framework that provides predictive performance estimates for many multilevel iterations and preconditioners. Early application of LFA was mainly focused on scalar problems or systems of PDEs with collocated discretizations. Standard LFA  cannot be directly applied to higher-order and staggered-grid discretizations, since the analysis depends on Toeplitz operator structure inherited from the mesh and discretization. There is a long history of generalization of “standard” LFA to account for more general structure of the discrete operator and/or multigrid algorithm, beginning with analysis of red-black (and, later, multicolour) relaxation [22, 33, 37]. This work has been generalized in recent years, leading to the idea of “periodic stencils”  and similar approaches for PDEs with random coefficients . From a different perspective, similar tools have arisen to account for the structure of coupled systems of PDEs, beginning with . This work was expanded by MacLachlan and Oosterlee , to account for the structure of overlapping relaxation schemes for the Stokes Equations. In a third setting, the mode analysis of certain space-time multigrid methods, Friedhoff and MacLachlan again use a similar approach to handle coarsening in time [18, 19]. In this work, we show how a generalization of these approaches can be applied to the structure of domain decomposition algorithms.
To our knowledge, there has been no research applying local Fourier analysis to BDDC-like algorithms. The same is true of the closely related finite element tearing and interconnecting dual-primal (FETI-DP) methodology [16, 17, 23]. Here, we adapt LFA to this domain decomposition method by borrowing tools from LFA for systems of PDEs and other block-structured settings. Noting that the stencil for the domain decomposition method depends on the size of the subdomains, an adaptation of the standard basis is useful, which we present in Section 4.1
. Because LFA can reflect both the distribution of eigenvalues and associated eigenvectors of a preconditioned operator, here, we adopt LFA to analyze variants of the common “lumped” and “Dirichlet” BDDC algorithms, based on, to guide construction of these methods. To do this, we use a modified basis as in [3, 21, 28] for the Fourier analysis that is well-suited for application to domain decomposition preconditioners.
Applying the two-level BDDC algorithm requires the solution of a Schur complement equation (coarse problem), which usually poses some difficulty with increasing problem size. Two- and three-level variants are, thus, considered in this paper. However, as is well-known in the literature, bounds on the performance of BDDC degrade sharply from two-level to three-level methods, particularly for large values of . Since our analysis shows that the largest eigenvalues of the preconditioned operator for the lumped BDDC algorithm are associated with oscillatory modes, we propose variants of standard BDDC based on multiplicative preconditioning and multigrid ideas. From the condition numbers offered by LFA, we can easily compare the efficiency of these variants. Furthermore, LFA can provide optimal parameters for these multiplicative methods, helping tune and understand sensitivity to the parameter choice.
This paper is organized as follows. In Section 2, we introduce the finite element discretization of the Laplace problem in two dimensions and the lumped and Dirichlet preconditioners. Two- and three-level preconditioned operators are developed in Section 3. In Section 4, we discuss the Fourier representation of the preconditioned operators. Section 5 reports LFA-predicted condition numbers of the BDDC variants considered here. Conclusions are presented in Section 6.
We consider the two-dimensional Laplace problem in weak form: Find such that
where is a bounded domain with Lipschitz boundary . Here, we consider the Ritz-Galerkin approximation over , the space of piecewise bilinear functions on a uniform rectangular mesh of . The corresponding linear system of equations is given as
We partition the domain, , into nonoverlapping subdomains, , where each subdomain is a union of shape regular elements and the nodes on the boundaries of neighboring subdomains match across the interface . The interface of subdomain is defined by . Here, we consider with both a discretization mesh (with meshsize ) and subdomain mesh (with meshsize ) given by uniform grids with square elements or subdomains.
The finite-element space can be rewritten as where is the product of the subdomain interior variable spaces . Functions in are supported in the subdomain and vanish on the subdomain interface . is the space of traces on of functions in . Then, we can write the subdomain problem with Neumann boundary conditions on as
is the restriction operator from a global vector to a subdomain vector on.
2.1 A Partially Subassembled Problem
In order to describe variants of the BDDC methods, we first introduce a partially subassembled problem, following , and the corresponding space of partially subassembled variables,
is spanned by the subdomain vertex nodal basis functions (the coarse degrees of freedom). The complementary space,, is the product of the subdomain spaces , which correspond to the subdomain interior and interface degrees of freedom and are spanned by the basis functions which vanish at the coarse-grid degrees of freedom. For a mesh, the degrees of freedom in are those corresponding to the circled nodes at the left of Figure 1, while the degrees of freedom in correspond to all interior nodes, plus duplicated (broken) degrees of freedom along subdomain boundaries.
The partially subassembled problem matrix, corresponding to the variables in the space , is obtained by assembling the subdomain matrices (3) only with respect to the coarse-level variables; that is,
where is a restriction from space to . The central idea behind BDDC preconditioners is to use as an approximation to that is more readily inverted, since the subdomain problems can be “condensed” to a solve on the underlying coarse grid, as described below.
The above “broken” space can, in fact, be viewed by considering the discretization of the PDE independently on each subdomain and only “assembling” the DOFs at the corners of the subdomains. Then, looking at this “subassembled” problem, the interface between any two subdomains has independent DOFs associated with each subdomain at the same locations, which are normally assembled to formulate the global system. This is different than the standard overlapping domain decomposition, which uses only one DOF at these locations in the overlap between subdomains and contributes this value to each of the shared subdomains. While can be readily defined in terms of , the reverse is not true.
2.2 Lumped and Dirichlet Preconditioners
In order to define the preconditioners under consideration for (2), we introduce a positive scaling factor, , for each node on the interface of subdomain Let be the set of indices of the subdomains that have on their boundaries. Define , where is the cardinality of . The scaled injection operator, , is defined so that each column of corresponds to a degree of freedom of the global problem (2). For subdomain interior and coarse-level variables, the corresponding column of has a single entry with value 1. Columns that correspond to an interface degree of freedom (the set of nodes in ) have non-zero entries each of .
Based on the partially subassembled problem, the first preconditioner introduced for solving (2) is
The preconditioned operator has the same eigenvalues as the preconditioned FETI-DP operator with a lumped preconditioner, except for some eigenvalues equal to 0 and 1 [17, 27]. We refer to as the lumped preconditioner. Note that corresponds to a naive approximation of by , using equal weighting of residual and corrections to/from the duplicated degrees of freedom in the partially broken space.
A similar preconditioner for
augments this by using discrete harmonic extensions in the restriction and interpolation operators, giving
where is the direct sum of local operators , which map the jump over a subdomain interface (given by ) to the interior of the subdomain by solving a local Dirichlet problem, and gives zero for other values. For any given , the component of on subdomain is given by
Extending the interface values using the discrete harmonic extension minimizes the energy norm of the resulting vector , giving a better stability bound. Furthermore, the preconditioned operator has the same eigenvalues as the BDDC operator , except for some eigenvalues equal to 1 . We refer to as the Dirichlet preconditioner. Note that differs from only by local operators in the operators used to map between the fully assembled and partially broken spaces.
Standard bounds (see, e.g., ) on the condition numbers of the preconditioned operators are that, for , there exists such that and, for , there exists such that .
While the preconditioners given by have some notational similarity to Galerkin corrections, it is important to note that acts on a higher-dimensional space than , and does not satisfy for either “restriction” operator. The essential difference between the two preconditioners is in how “glues” the solution over the broken space back into the usual continuous space. does this through a simple partition of unity, while minimizes subdomain energy when propagating mismatches along the subdomain boundaries.
3 Two- and Three-level Variants
In both of the above preconditioned operators, we need to solve the partially subassembled problems, that we write in block form
is the identity matrix,and contains the subdomain interior and interface degrees of freedom, and corresponds to the coarse-level degrees of freedom, which are located at the corners of the subdomains. We write in (8) in factorization form to easily separate the action on the coarse degrees of freedom, and to find the corresponding symbol of . If we define
then the Schur complement is simply the Galerkin coarse operator, , and the block-factorization solve for is equivalent to a two-level additive multigrid method with exact -relaxation using
to define an “ideal” relaxation scheme to complement the coarse-grid correction defined by .
In the partially subassembled problem (8), we need to solve a coarse problem related to . We can either solve this coarse problem exactly (corresponding to a two-level method, where the Schur complement is inverted exactly) or inexactly (as a three-level method), where the lumped and Dirichlet preconditioners defined above are used recursively to solve this problem.
3.1 Exact and Inexact Solve for the Schur Complement
and note that and are formed from the LDU factorization of as in (8). For , let denote the preconditioned operators for two- and three-level variants of BDDC, where and denote using and (with as preconditioners for the fine and coarse problems, respectively, where stands for applying the preconditioner to the Schur complement matrix, . By standard calculation, we can write
When , is a two-level method, solving the Schur complement problem exactly, as . Note that, for the three-level variants ,
corresponding to replacing in (8) by representing the inexact solve. From this, it is clear that can be applied without directly applying the inverse of .
The eigenvalues of are real and bounded below by 1.
which is symmetric and positive definite (SPD), ensuring the eigenvalues of are real. Similarly, can be rewritten as
Writing , we have
which leads to
Now, for any , taking in (9), we have
According to [27, Theorem 4], we then have
which means that the smallest eigenvalue of is not less than 1.
3.2 Multiplicative Preconditioners
As we shall see, the bounds above are relatively sharp and the performance of both preconditioners degrades with subdomain size and number of levels. To attempt to counteract this, we consider multiplicative combinations of these preconditioners with a simple diagonal scaling operator, mimicking the use of weighted Jacobi relaxation in classical multigrid methods. We use to denote the multiplicative preconditioned operator based on with diagonal scaling on the fine level. Here,
where is the diagonal of and is a chosen relaxation parameter. Note that , so represents the multiplicative combination given.
The eigenvalues of are real.
First, recall that is SPD, so it has a unique SPD matrix square root, . Note that . From Theorem 1, we know that is symmetric, so is symmetric as well.
is similar to , which can be rewritten as
From Theorem 1, we know that the eigenvalues of are not less than 1. Thus, is symmetric positive semi-definite, and has a unique symmetric positive semi-definite square root. Using the fact that , we have
Note that is symmetric. Thus, the eigenvalues of are real and so are those of .
Another variant is the use of multiplicative preconditioning on the coarse level with a similar diagonal scaling. We use to denote the resulting multiplicative preconditioner. Here,
where is the diagonal of .
Instead of using a single sweep of Jacobi in , we can consider a symmetrized Jacobi operator , where ; that is,
then changes to
for defined as but with in the -block.
Finally, we can also apply the multiplicative operators based on diagonal scaling on both the fine and coarse levels. We denote this as
where is the diagonal of and is a chosen relaxation parameter.
The eigenvalues of operators , and are not generally all real. When , is similar to a symmetric matrix and we observe a real spectrum in numerical results, but not for all cases when . In most situations, and have complex eigenvalues.
In the following, we focus on analyzing the spectral properties of the above preconditioned operators by local Fourier analysis . The main focus of this work is on the operators and , because the Fourier representations of other operators are just combinations of these three and some simple additional terms.
4 Local Fourier Analysis
To apply LFA to the BDDC-like methods proposed here, we first review some terminology of classical LFA. We consider a two-dimensional infinite uniform grid, , with
and Fourier functions on , where and . Let be a Toeplitz operator acting on as
with constant coefficients , where is a function on . Here, is taken to be a finite index set. Note that since is Toeplitz, it is diagonalized by the Fourier modes .
If for all grid functions, ,
we call the symbol of .
In Definition 1, the operator acts on a single function on , so is a scalar. For an operator mapping vectors on to vectors on , the symbol will be extended to be a matrix.
4.1 Change of Fourier Basis
Here, we discuss LFA for a domain decomposition method. While the classical basis set for LFA, denoted below, could be used, we find it is substantially more convenient to make use of a transformed “sparse” basis, introduced here as . This basis allows a natural expression of the periodic structures in domain decomposition preconditioners that vary with subdomain size. We treat each subdomain problem as one macroelement patch, and each subdomain block in the global problem is diagonalized by a coupled set of Fourier modes introduced in the following. Because each subdomain has the same size, , we consider the high and low frequencies for coarsening by factor , given by
Let , where and for . For any given , we define the -dimensional space
as the classical space of Fourier harmonics for factor coarsening.
For any , we consider a grid function defined as a linear combination of the basis functions for with frequencies and coefficients as
We note that any index has a unique representation as where and . From (15), we have
Thus, we can write
In other words, for any point with and , can be reconstructed from a single Fourier mode with coefficient . Thus, on the mesh defined in (14), the periodicity of the basis functions in can also be represented by a pointwise basis on each -block.
Based on (16), we consider a “sparse” -dimensional space as follows
Note that, with this notation, (16) can be rewritten as
and are equivalent.
While the derivation above shows directly that , we revisit this calculation now to show that the mapping is invertible and, hence, as well.
Let be an arbitrary vector with size , denoted as
Then, we define a vector, , based on (17), as follows
Let be the matrix of this transformation, , and
Note that defines a vector whose -th entry is and, thus, we see that .
Note that is a Vandermonde matrix based on values , where . It is obvious that if . Consequently, . Thus, is invertible, and so is . It follows that and are equivalent.
Let , be the primitive -th root of unity, and note that . Thus,
is the unitary discrete Fourier transform (DFT) matrix with