{'fraction_non_alphanumeric': 0.035016510731975786, 'fraction_numerical': 0.012887543569987158, 'mean_word_length': 4.493009195112735, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A practical factorization of a Schur complement for PDE-constrained Distributed Optimal Control

A distributed optimal control problem with the constraint of a linear elliptic partial differential equation is considered. A necessary optimality condition for this problem forms a saddle point system, the efficient and accurate solution of which is crucial. A new factorization of the Schur complement for such a system is proposed and its characteristics discussed. The factorization introduces two complex factors that are complex conjugate to each other. The proposed solution methodology involves the application of a parallel linear domain decomposition solver---FETI-DPH---for the solution of the subproblems with the complex factors. Numerical properties of FETI-DPH in this context are demonstrated, including numerical and parallel scalability and regularization dependence. The new factorization can be used to solve Schur complement systems arising in both range-space and full-space formulations. In both cases, numerical results indicate that the complex factorization is promising.

1. Introduction. Numerical methods for solving partial differential equations (PDEs) have broad applications in the simulation of complicated physical models, the prediction of their response, and design. An important and practical subset of these applications-particularly for design-involves the use of a mathematical optimization technique in which the PDE takes the role of a constraint equation. Two methodologies for PDE-constrained optimization problems are SAND (simultaneous analysis and design) and NAND (nested analysis and design) [3]. NAND uses PDEs to express decision variables as an implicit function of state variables and does not include state variables as optimization variables. Thus, the size of the optimization problem is not typically large. On the other hand, the SAND approach takes both decision variables and state variables as optimization variables and considers PDEs to be equality constraints. Consequently, the size of the system in the SAND approach is generally much larger. The NAND approach has traditionally been the method of choice for physics-based applications not only because SAND requires the solution of a large-scale system of equations but also because NAND conveniently permits the direct use of existing solvers for both optimization and PDE simulation as a black box. However, NAND suffers from the fact that many PDE simulations are required for function evaluations-typically the most expensive part of this approach. Due to the continuous increase in computational power (e.g., speed of processing, capacity of memory, high performance computing) accompanied by the development of robust and versatile numerical algorithms (e.g., parallel algorithms), the SAND approach has received increasing attention from researchers in recent years [4,5,7,31], and the present study continues this line of research.
The SAND approach to PDE-constrained optimization takes the form minimize y,u F (y, u) subject to C(y, u) = 0, (1.1) where C(y, u) = 0 is the time-independent PDE constraint, y is the vector of state variables, and u is the vector of decision variables. State variables by definition are the unknown variables in the forward PDE problem. For example, state variables comprise temperature in heat conduction problems and displacements in elastostatics. For the class of PDE-constrained optimization known as optimal control the decision variables u are referred to as control variables, whereas for optimal design or shape optimization problems the decision variables u are called design variables. The decision variables may also be a set of parameters describing the material properties or the system of dynamics for some inverse problems. All three types of PDE-constrained optimization problems share a similar structure of their linear or linearized system of equations known as a KKT system, after the Karush-Kuhn-Tucker optimality conditions, or a saddle point system.
In what follows, a robust and versatile numerical method for solving a distributed optimal control problem is considered. In particular, heat conduction and elastostatic PDE-constrained problems are studied, although the method developed here can also be applied or extended easily to other types of PDE-constrained optimization problems. From various possible objective functions that may be used to formulate a distributed PDE-constrained optimal control problem, the one considered is that in which a target state is assumed to be given. Thus, the aim is to find a state that is close to the prescribed target and a control that realizes that particular state. For example, the static thermal conduction optimal control problem with a target temperature distributionȳ is formulated as minimize y,u F (y, u) := 1 2 Ω (y −ȳ) 2 dx + φ 2 Ω u 2 dx subject to −∇ 2 y = u on Ω y = y c on Γ g . (1. 2) The solution of this problem is considered in Section 6.1. Variables y and u are the temperature state and heat source control, respectively, while y c is the prescribed boundary conditions and φ is a regularization parameter. The domain of interest is denoted as Ω and a Dirichlet boundary condition is imposed on the boundary Γ g . Note that a unit conductivity is assumed. Two alternative approaches to PDE-constrained optimization problems such as (1.2) are optimize-and-discretize and discretize-and-optimize. The latter approach, in which one first discretizes both the objective function and constraints and then obtains the discretized optimality condition, is typically used for PDE-constrained optimal control [5,9,30,29,32,35,11] and is followed here.
Solving the saddle point system representing the discretized optimality condition efficiently is crucial to the competitiveness of the SAND approach, in comparison with NAND. There are two methodologies for solving a saddle point system.
• In reduced-space methods, one attempts to reduce the size of a saddle point system by eliminating some variables and solving a smaller system [37,38,35,4]. The range-space method and the null-space method are two popular reduced-space methods. In the range-space method, one solves for the dual variables first using the corresponding Schur complement and subsequently updates the primal variables, whereas the null-space method subdivides the variables algebraically into null-space variables and range-space variables using null-space and range-space bases of the constraint Jacobian. • In full-space methods, one solves for the primal and dual variables simultaneously. The resultant system of equations becomes a sparse saddle point system and iterative methods are the only practical choice for large problems. Various efficient preconditioners for saddle point systems have been developed [27,2,22,16,32,30,10,9,23,5,36]. However, there is still motivation for further research in this area. For example, most preconditioners are not robust when φ is small. Recently, Pearson and Wathen have developed a new approximation of the Schur complement and used it to facilitate a regularization-robust preconditioner for a particular distributed optimal control problem [30]. Pearson et al. have further extended the usage of the aforementioned approximation to a broader range of optimal control problems [29].
The factorization of the Schur complement presented in this paper has a similar form to Pearson and Wathen's approximation, and furthermore is applied to the same particular distributed optimal control problem [30]. However the approach taken here differs from that of Pearson and Wathen in two ways. First and foremost, because the factorization is exact by its nature, the range-space method can be adopted and one efficient solve with the Schur complement results in a solution to the distributed optimal control problem. Second, a scalable domain decomposition based parallel linear solver FETI-DPH [13] is used to solve each of the subproblems arising in the application of the proposed factorization. In contrast, Pearson and Wathen suggest a multigrid method.
The outline of the paper is as follows: Section 2 describes discretization of the distributed optimal control problem (1.2) and the corresponding optimality condition, which is a saddle point system. Section 3 outlines two methods for solving the saddle point system: the range-space method and the full-space method. Additionally, this section reviews existing preconditioners for the full-space method related to the Schur complement in the range-space method. Section 4 introduces a practical factorization of the Schur complement and explains how it can be used in both the range-space method and the full-space method. Section 5 describes the Finite Element Tearing and Interconnecting (FETI) method and its so-called dual-primal variants, FETI-DP and FETI-DPH. Section 6 presents numerical results that illustrate the scalability and efficiency of the method applied to a selection of distributed optimal control problems.
2. Distributed control problems. In this section, the finite element discretization of a distributed optimal control problem is introduced and the corresponding optimality condition is presented. For brevity, we choose to show the discretization of a distributed optimal control problem with a constraint of the Poisson equation (i.e., Eq. (1.2)). However, distributed optimal control problems with other types of PDE constraints (e.g., linear elasticity) will have the same discretized formulation. The finite element discretization of (1.2) gives minimize y,u where K, V ∈ R n×n , and K c ∈ R n×m are the stiffness matrix, volume matrix, and constrained stiffness matrix, respectively. Vector valued quantitiesȳ, y, u ∈ R n , and y c ∈ R m are the discretized versions of the target state, state, control, and prescribed boundary conditions, respectively. All the discrete variables are denoted with bold fonts for the remainder of the paper. The dimensions n and m are the number of unconstrained degrees of freedom (i.e., the size of y) and the number of constrained degrees of freedom, respectively. It is assumed that both K and V are symmetric positive definite (SPD) matrices, which is the case for the thermal problem above. The energy norm · V is defined as q V = q T Vq. The details of the finite element discretization procedure can be found in [7]. Problem (2.1) is a convex quadratic program and its solution is a saddle point of the Lagrangian, A saddle point (y * , u * , λ * ) must satisfy the following optimality condition: This equation is of the form Ax = b, where A is a symmetric indefinite matrix. Section 3 explains two ways of solving the saddle point system in (2.3): the rangespace method and full-space method.
3. Range-space and full-space methods.
3.1. Range-space method. The range-space method was introduced for largescale inequality-constrained convex quadratic programming problems with small active sets in order to overcome the disadvantage of the null-space method: the increasingly high dimension of the null space basis as a solution is approached [19,20]. For equality-constrained convex quadratic programming, the range-space method is equivalent to the Schur complement method [37] where the corresponding dual problem is solved first. The range-space method is a reduced-space method in a sense that the size of the saddle point system in (2.3) is reduced by eliminating y * and u * and solving for λ * first. The dual variable λ * is obtained by solving where S = KV −1 K + 1 φ V is known as the negative Schur complement on the dual variables. Then, y * and u * are computed from The most expensive and crucial step in the range-space method is to solve with S in (3.1). The eigenvalues (and consequently, the condition number) of S depend on those of K and V and the value of φ [37,7]. For small values of φ, the condition number of S is affected most by the eigenvalues of V, whose condition number is bounded above by a constant, C, for P q or Q q (the q th order triangular or quadrilateral finite elements in 2D and tetrahedral or brick elements in 3D) if a set of grids is quasiuniform (see Eq. (1.116) and Eq. (1.117) in [12]). The constant C depends on the order of approximation q but not on the mesh size h. Thus for small values of φ, even if the grid is refined or the problem size is increased, the condition number of S is bounded. For large values of φ, the condition number of S (denoted by κ(S)) depends on κ(KV −1 K), which is proportional to κ(K) 2 . For a second-order elliptic operator, one can prove that κ(K) < ch −2 if P q or Q q is used on a quasi-uniform discretization of a domain (see Eq. (1.119) and Eq. (1.121) in [12]). For a higher-order elliptic operator, the dependence on h of κ(K) increases (e.g., the shell or beam elements). This is unfortunate because κ(K) is likely to increase as the mesh size h decreases or the problem size increases. One remedy is to apply iterative refinement with a high precision data type. If the solution is inaccurate even after iterative refinement due to the ill-conditioning of S, an alternative way of solving (2.3) that avoids an exact solve with S is to use the full-space method, which is explained in the next section. An additional challenge in applying the range-space method for solving (2.3) resides in solving with S itself. Because S is a sum of two matrices (KV −1 K and 1 φ V), the first of which is a product requiring two matrix-matrix multiplications to evaluate, it is not straightforward to come up with an efficient way of solving with S. In Section 4, a practical factorization of S is introduced and the factorization does not require matrix-matrix products. The factorization introduced in Section 4 introduces complex symmetric matrices (not Hermitian matrices) in its factors. A list of efficient solvers for complex symmetric matrices includes CG-type methods [17], CS-MINRES-QLP [6], GMRES [34,8], FETI-DPH [13], and multi-grid methods [24,33]. FETI-DPH is explained in Section 5 and used in numerical experiments in Section 6.

3.2.
Full-space Method. The full-space method attempts to solve for all the variables of the saddle point system in (2.3) simultaneously. As the discretization is refined, the size of the system becomes large and an iterative method is often the only available method. Because A is a symmetric indefinite matrix, any Krylov iterative method suitable for this class of matrices, such as MINRES [28], SYMMLQ [28], SQMR [18], and GMRES, can be used. For successful convergence of the iterative method, one is often required to apply a good preconditioner. Even without any preconditioner, the saddle point system itself tends to be better conditioned for moderately large values of φ than the Schur complement in the range-space method [37]. If a good preconditioner is used, the full-space method is likely to have a better scalability than the range-space method, meaning that computational cost does not grow at an exponential rate as the problem size increases. Thus, many preconditioners have been developed recently. Here we focus on a Schur complement-based preconditioner because the primary interest is to discuss the usage of a Schur complement factorization that will be introduced in Section 4. Murphy, et al. [27] have shown that if the following block diagonal preconditioner is used, then the preconditioned system has at most three nonzero distinct eigenvalues: where S is the negative Schur complement defined in (3.1). This implies that the maximum number of iterations required for convergence in a Krylov iterative method is three if P mgw is used as a preconditioner. However, applying this preconditioner has previously been considered impractical because it requires solving with S in order to apply P mgw . Even if it were practical to solve with S efficiently and accurately, then the preferred approach would typically be to use the range-space method rather than to apply P mgw in the full-space method because the range-space method requires only one solve with S, while the full space method with P mgw as a preconditioner is likely to require more than one solve with S. There has been some research done on finding a good approximation of the Schur complement for use in a Schur complement-based preconditioner (e.g., P mgw ) [32,30,29]. Particularly, Pearson and Wathen [30] have developed the following approximation, which is regularization-robust: 4) and proved that the eigenvalues of S −1 p S are between 1 2 and 1 regardless of the mesh size h and φ. This implies theoretically that a Krylov iterative method must converge in O(1) iterations regardless of h and φ. They have used an algebraic multigrid method to solve with the factor K + 1 √ φ V in their numerical examples and demonstrated that the number of iterations required for convergence is indeed O(1). However, in the problems they solve, the condition number of the Schur complement is small enough so that the approximation to the Schur complement behaves well. The factorization introduced in Section 4 has a similar form to (3.4). However, it is an exact representation of S and therefore enables the application of the range-space method.

4.
A "practical" factorization of the Schur complement. The Schur complement S in (3.1) can be factored into the following form: The form of this factorization is similar to S p in (3.4) in the sense that the first and third factors are some linear combinations of K and V, and the middle factor is V −1 . This suggests that the same methods for solving with the factors in S p may also be used to solve with factors in S c . However, there are some important differences between S p and S c . The first and third factors of S p are the same, but the corresponding factors of S c are not. Consequently, if direct solvers were to be used for each factor, then only one factorization would be required for S p , while two factorizations of two different systems would be required for S c . Additionally, S c introduces complex numbers, but all the elements in S p are real. Thus, one solve with S c requires four times more storage and floating point operations than one solve with S p . In spite of the disadvantages of applying S c , it is an exact representation of S, and this permits use of the range-space method where only one solve with S c is sufficient to obtain a solution to (2.1). On the other hand, many solves with S p are required because S p is an approximation to S, so it can only be used as a preconditioner. This discussion leads to the following two extreme cases when dealing with the Schur complement: • If a direct method is the preferred choice for solving the Schur complement (i.e., a factorization of a given system is required), it would be advantageous to apply the direct method to S p assuming that the dominating computational cost occurs in factorization of the system. • If an iterative method is the only option and the range-space method is applicable, then it would be preferable to apply the iterative method to S c . Between these two extreme cases, a choice must be made depending on the characteristics of the problem and a numerical solver. For example, if the domain of a problem is complex, then the computational overhead incurred by applying S c rather than S p will be substantially diminished since in this case the factors of both S c and S p will be complex. Such problems include any frequency domain analyses with damping in structural or acoustic problems. On the other hand, if a problem requires many solves with a Schur complement and with multiple right-hand sides and is small enough to permit a direct solver to be used for the factors of S, applying S p is favorable because direct methods are more efficient than iterative methods in general for multiple righthand sides. If a domain decomposition method-in which both direct and iterative methods are used-is chosen, then one needs to examine the costs of each component of the computation (e.g., data structure, building and storing the operator, and solving with the operator). In Section 6, numerical results for solving with S c are shown using the domain decomposition solver FETI-DPH [13].

FETI.
In order to be self-contained, the FETI method [15] and two of its variants (FETI-DP [14] and FETI-DPH [13]) are explained in this section. For a more detailed description, see [15,14,13,1]. The FETI method was developed in order to solve, in parallel, the following linear system of equations arising from the finite element discretization of a linear elasticity PDE: where K is the stiffness matrix (or a linear combination of the stiffness matrix and mass matrix), y is the vector of unknown displacements, and f is an external force term. Solving (5.1) is equivalent to minimizing a quadratic function: The objective function in (5.2) often represents a physical quantity (e.g., energy in linear elasticity). In domain decomposition methods, the spatial domain Ω is divided into N s either overlapping or non-overlapping subdomains. The FETI method divides Ω into N s non-overlapping subdomains as illustrated in Figure 5.1 and minimizes a local function in each subdomain (i.e., minimize y sT K s y s − y sT f s in Ω s where the superscript "s" designate the restrictions to the specific subdomain). In order to be equivalent to (5.2), the following additional continuity condition on y s between interfaces must be imposed: where B s is a signed Boolean matrix. Hence one can write the following constrained quadratic programming equivalent to (5.  Statement (5.6) says that one can find an optimal solution y s * in two ways. The first way is to minimize L(y s , λ) with respect to y s , then maximize the resultant with respect to λ. The second way is to maximize L(y s , λ) with respect to λ, then minimize the resultant with respect to y s . The FETI method takes the first approach where the Lagrange dual function g(λ) = inf y s L(y s , λ) is obtained and then maximized in order to obtain the Lagrange multipliers first. This approach is attractive if the number of constraints (i.e., interface continuity condition) is small, which is expected to be the case in general.
For a fixed λ, L(y s , λ) is separable in y s (i.e., 1 2 y sT K s y s − y sT f s + λ T B s y s on Ω s ). Thus, the Lagrange dual function can be obtained by solving K s y s = f s −B s λ on each subdomain Ω s . In the case of K being a stiffness matrix as in (5.1), K s is symmetric positive definite or semidefinite. If K s is positive semidefinite, then it is necessary to explicitly ensure that K s y s = f s − B s λ is compatible (i.e., R sT (f s − B sT λ) = 0 where the columns of R s span the left null space of K s ). Otherwise, the minimum of 1 2 y sT K s y s − y sT f s + λ T B s y s in terms of y s becomes −∞ and the objective value of primal problem (5.4) is also −∞ by strong duality. However, this is not a physical solution. Thus, we restrict ourselves to the case when g(λ) > −∞.
Finally, the Lagrange dual function g(λ) = inf y s L(y s , λ) is defined as and K s+ is a pseudo-inverse of K s . This defines the following dual problem to (5.4): The FETI method solves (5.9) with a Preconditioned Conjugate Projected Gradient (PCPG) algorithm. The FETI formulation above is suitable for parallel processing. Any subdomain level computations (e.g., factorization or computation with K s and f s ) can be assigned to an individual process. The size of λ is the total number of degrees of freedom restricted to the interfaces between subdomains. Within PCPG, λ is projected onto the domain of feasibility, which in turn reduces the effective dimension of the problem. Indeed the FETI method equipped with the Dirichlet preconditioner is proven numerically scalable with respect to both problem size and number of subdomains. The Dirichlet preconditioner is defined as 11) and the subscript i and b denote subdomain internal degrees of freedom and interface degrees of freedom, respectively. If it is applied in the conjugate gradient algorithm with the projected gradient [21] to second-order elliptic problems, the condition number κ of the interface problem (5.9) is approximately [26] However, the first-generation FETI method is not numerically scalable for fourthorder plate and shell problems. This leads to FETI-DP. FETI-DP is one variant considered to be the third-generation FETI method. FETI-DP enforces continuity at some interface corner nodes at each iteration (see Figure 5.1). An extra coarse problem needs to be solved [14] and in each subdomain the remaining subdomain stiffness matrix K s rr (i.e., that excluding the degrees of freedom in the corner nodes at which the continuity is enforced) becomes positive definite. Consequently, the corresponding dual interface problem becomes an unconstrained QP. This treatment helps to achieve numerical scalability for fourth-order elliptic problems. The continuity constraints can be augmented by additional constraints that are enforced exactly throughout the iterations in FETI-DP in order to accelerate the convergence. This augmentation procedure results in the augmented coarse problem in FETI-DP [25]. A standard augmentation procedure uses the edge-based rigid body modes (rotational and/or translational) [13]. The positive definiteness of K s rr is not guaranteed when K in (5.1) is different from the stiffness matrix. For example, for a Helmholtz problem or frequency response elastodynamic problem, K in (5.1) becomes where C is a symmetric matrix that arises from the discretization of an absorbing boundary condition, and k > 0 is a frequency (or a wave number for acoustic scattering problems). In this case, depending on the value of k, Z may become indefinite. The same difficulty is encountered when solving with complex factors in (4.1). In order to resolve this difficulty, two special treatments were required.
• The first treatment is required to deal with indefinite matrices. A solver that is suitable for indefinite matrices must be applied. Such an iterative solver includes GMRES [34] for a general square matrix and MINRES [28] for a symmetric indefinite matrix. FETI-DPH ('H' stands for Helmholtz) [13] uses GMRES to deal with indefiniteness. • The second treatment is required for improved convergence when the wave number of frequency is large. A particular augmentation used in FETI-DPH that accelerates convergence for the Helmholtz or frequency response elastodynamic problem is the free-space solutions of the corresponding equations, which are plane waves. Note that the complex factors in Section 4.1 are similar to (5.13), but not identical, and not surprisingly the plane wave augmentation was found to not provide the same beneficial effects as in the Helmholtz problem. Thus, the edge-based rigid body modes (both rotational and translational) augmentation is applied in the next section for two numerical examples. Strictly speaking, although FETI-DPH normally refers to the variant of FETI-DP that uses both GMRES and plane wave augmentation, in what follows, FETI-DPH refers to a FETI-DP solver incorporating GMRES and edge-based rigid body modes.
6. Numerical Results. Two pedagogical examples are considered: linear static heat control of a 2D square plate and structure control of a 3D solid cantilever. The heat control problem in Section 6.1 is taken identical to the thermal problem solved in the paper by Rees et al. [32]. The Schur complement that arises in this thermal problem is sufficiently well-conditioned for the range of φ values and mesh size h considered here for the range-space method to be applicable and therefore it is the only method used. The cantilever control problem in Section 6.2 is very flexible with material properties of rubber. For a relatively large φ value, the range-space method fails to converge to an accurate solution. Thus, the range-space method is used for small values of φ only, and the full-space method is used for large values of φ. In the case of the full-space method, both the approximate Schur complement proposed by Pearson and Wathen [30] and the complex factorization introduced in this paper are used in a Schur-complement based preconditioner and compared in terms of computational time and iteration counts.
It is well established that augmentation is a beneficial feature that improves the performance of FETI. However, the optimality of the augmentation's performance depends on the nature of the problem. Because neither the FETI solver nor any of its variants has ever been applied to a system of the form K ± 1 i √ φ V, the numerical performance of FETI-DPH with and without augmentation is compared in this section. All the simulations were run on a heterogeneous Linux cluster with 2.6 GHz hexa-core Westmere processors (32 blades, 2 processors per blade, QDR Infiniband interconnect, 24GB/blade) and 2.6 GHz octa-core Sandybridge processors (4 blades, 2 processors per blade, FDR Infiniband interconnect, 256 GB/node).

Thermal problem.
A linear static heat control problem is solved. This example is identical to example 5.1 in [32] by Rees et al. The same objective function and PDE constraints are used as in (1.2), whose continuous formulation is reproduced here for better access: The domain Ω is [0, 1] 2 ⊂ R 2 , which is a unit square plate, whose heat conductivity is unity. The target temperatureȳ is defined as which is illustrated in Figure 6.1. The boundary condition y c is defined as The optimal control problem (6.1) tries to find a temperature y that is close to the target temperatureȳ by controlling heat u. How close y can be toȳ is determined by the parameter φ. As φ decreases, y is expected to approachȳ, but u is also expected to increase. This makes sense because the objective function value is not sensitive to the second term if the parameter φ is small and the first term dominates the objective function value. The results shown in Figure 6.2 confirm the validity of these expectations. For this particular example, according to Figure 6.2(a), one needs φ less than 2 × 10 −6 in order to match the target temperature to within 1%.  According to Figure 6.2(b), however, one needs to set φ greater than 2 × 10 −4 in order for the norm of control solution not to exceed 100.
The same effects can be demonstrated visually in Figure 6.3. Temperature and heat distributions for various φ values are shown in Figures 6.3(a)-(c). Note that for φ = 0.02, the heat is almost zero everywhere and the corresponding temperature distribution (induced mainly by the boundary condition y c ) is slightly different from the target temperature (i.e., Figure 6.1). The target temperature can be matched more precisely if a smaller φ is used. Figures 6.3(b) and (c) show temperature distributions that are closer to the target temperature. They are produced by using smaller φ values (i.e., φ = 2 × 10 −4 and 2 × 10 −5 ). Note that the corresponding heat distributions show noticeable non-zero heat at the left bottom of the domain. As φ decreases, a sharper heat gradient is visible near the boundary. Table 6.1 shows FETI-DPH's dependence on φ and the effects of augmentation   Table 6.1 shows that, in the absence of augmentation, the number of iterations tends to increase, but not strictly, as φ decreases. This dependence on φ is alleviated by the introduction of rigid body modes augmentation and some CPU time is saved. Table 6.1 demonstrates that rigid body modes augmentation works well for the thermal optimal control problem. The effects of rigid body modes augmentation is even more dramatic in the case of solid elements, as described in the following section. Further research on the existence of an optimal augmentation for the thermal optimal control problem is an interesting future topic (e.g., a set of modes that decreases the number of FETI-DPH iterations to an order of a constant regardless of the values of φ). Table 6.1 also shows that almost the same number of iterations of FETI-DPH is required on This result is not surprising, given that the two operators are complex conjugate to each other, and consequently have complex conjugate eigenvalues, the same pattern of clusterings, and an identical condition number.
There are two kinds of scalability tests for domain decomposition methods: numerical scalability and parallel scalability [13]. Numerical scalability is studied by varying the problem size, subdomain size, and the number of subdomains. Table 6.2 shows numerical scalability of FETI-DPH on the thermal problem (6.1). The number of iterations is shown for various mesh sizes h and for various subdomain sizes H. The regularization parameter φ is 2 × 10 −8 . The convergence threshold for FETI-DPH is 10 −10 . The fixed ratio H/h = 2 6 is used. Because the ratio H/h is fixed to be 2 6 , each subdomain contains 4096 elements. The problem size is varied from around 15 thousand degrees of freedom to 67 million degrees of freedom. Based on the theoretical result of (5.12), the number of FETI-DPH iterations must be more or less constant no matter what problem size is considered, provided that the ratio H/h is fixed. Indeed, Table 6.2 shows that it is the case. As h decreases, the total number of elements increases, meaning that the problem size increases. Table 6.2 shows that the number of FETI-DPH iterations actually decreases as the problem size increases, which is consistent with (5.12). Again, almost the same number of iterations of FETI-DPH is Parallel scalability measures how fast a domain decomposition method converges for a fixed problem size, a fixed subdomain size, and a fixed number of subdomains with increasing number of processes. Table 6.3 presents the parallel scalability of FETI-DPH for the heat conduction control problem (6.1). A mesh size h = 1/1024 (around 1 million degrees of freedom) and subdomain size H = 1/16 (4096 elements in each subdomain) are used. The regularization parameter φ is 2 × 10 −8 . The convergence threshold is 10 −10 . The number of processes (N p ) varies from 1 to 144 and parallel speed-up is measured relative to the FETI-DPH CPU time of N P = 1. As the number of processes increases up to 144, Table 6.3 shows that a parallel speed-up of 26.8 is gained. The thermal element (i.e., finite element for Laplacian operator) is widely used when one wants to verify a new numerical method or to analyze it. In the next section, a more difficult problem is considered (in the sense that properties of a real material are used and 3D solid elements are used).

Solid Cantilever Control.
In this section, the proposed optimal control techniques are applied to a linear static structural problem, a cantilever with solid elements. In this case, the constraints in (6.1) are replaced by an elastostatic PDE. The cantilever has a square cross-section of 1 × 1 m 2 and a length of 3 m. The target displacementȳ is generated by running a forward PDE simulation with a uniform pressure load of 100 kPa applied to the bottom surface of the cantilever. The fol- . Note that the left end is completely fixed. The body force is used as a control variable. Although this control variable is not practical because the body force control is not physically attainable, this problem is useful to investigate how the FETI-DPH solver performs on each factor of the Schur complement factorization for problems involving solid finite elements and linear elastostatic PDEs.   i √ φ V vary as φ varies. Both condition numbers increase as φ increases but the order of magnitude of the condition number of the factor is one half that of S. As φ increases, the accuracy of the range-space method degenerates. For example, for values of φ larger than 10 −2 , the constraint violation becomes larger than 10 −4 , which is of questionable acceptability. On the other hand, the full-space method shows high accuracy consistently for the entire range of φ values considered. It is possible to achieve small improvements in the accuracy of the range-space method by using iterative refinement and higher floating-point precision. However, this approach incurs an additional cost and preliminary experiments suggest that it is not a competitive solution. One could argue that the values of interest of the parameter φ for this problem are the ones smaller than 10 −5 because the relative difference between target and solution states is greater than 83% for values of φ outside this range. For this range of φ values, the accuracy of the range-space method is acceptable for this particular mesh size (h = 1/8).  Table 6.4 shows FETI-DPH's dependence on the value of φ for the cantilever control problem. The mesh size is h = 1/90 and the convergence threshold is 10 −9 . The size of the subdomain is H = 1/6 and the number of processes is 64. The second and third columns show the number of FETI-DPH iterations without any augmentation, and the fourth column shows the corresponding FETI-DPH CPU time in seconds. The signs "-" and "+" indicate the FETI-DPH solves on K − 1 i √ φ V and K + 1 i √ φ V, respectively. The fifth and sixth columns show the number of FETI-DPH iterations with edge-based rigid body modes augmentation, and the seventh column shows the corresponding FETI-DPH CPU time in seconds. The last column shows the speedup due to augmentation. Without augmentation, the performance of the FETI-DPH solver degrades as the value of φ increases from 10 −12 to 10 −7 . This makes sense because the condition number of the complex factors in (4.1) and the Schur complement itself increase as φ increases (see Figure 6.5). However, for the range of smaller φ values (i.e., less than 10 −12 ), the number of FETI-DPH iterations without augmentation increases as φ decreases. This is analogous to FETI-DPH being dependent on wave numbers in acoustic problems if no augmentation is applied (i.e., more iterations are required for a larger wave number without augmentation) although the complex factors are not the same as in the system that normally arises from a Helmholtz problem. In order to alleviate this deterioration, edge-based rigid body mode augmentation is used. Table 6.4 shows that the augmentation decreases the number of FETI-DPH iterations substantially and reduces the CPU time by a factor of four on average. However, the rigid body modes augmentation is not optimal in the sense that as φ decreases the number of FETIP-DPH iterations still increases. Table 6.5 shows the numerical scalability of the FETI-DPH solver for the cantilever control problem. The number of iterations is shown for various mesh sizes h and for various subdomain sizes H. The regularization parameter φ is 1 × 10 −16 . The convergence threshold for FETI-DPH is 10 −9 . The problem size is varied from around 2 million degrees of freedom to 40 million degrees of freedom. As in the numerical scalability test for the thermal problem, the H/h ratio is fixed to 15 (exactly 3375 elements in each subdomain) and the number of FETI-DPH iterations is counted for As the problem size increases, the number of FETI iterations required for convergence slightly decreases, which is again consistent with the theoretical result (5.12).  Table 6.6 shows the parallel scalability of FETI-DPH for the cantilever control problem. The FETI-DPH CPU time and the corresponding parallel speed-up are shown for the various numbers of processes. The regularization parameter φ is 10 −12 . The convergence threshold is 10 −9 . A mesh size h = 1/90 (i.e., around 6.6 million degrees of freedom) and a subdomain size H = 1/6 are used. The number of processes (N p ) increases from 1 to 64 and parallel speed-up is measured as the ratio with respect to the CPU time of N P = 1. As the number of processes increases up to 64, Table 6.6 shows the speed up is greater than half the maximum possible.  (3.4)). Thus, one must consider the effectiveness of the complex factorization carefully although only one solve with complex factorization is required in the range-space method and fewer iterations are required in the full-space method. Table 6.7 reports comparison of CPU time and number of iterations between S p and S c in the full-space method. The block diagonal preconditioner of Murphy, et al. (3.3) is used. Comparison is made for various values of φ and each solve is done by FETI-DPH with edge augmentation. The mesh size is h = 1/30 and the convergence threshold for FETI-DPH is 10 −9 . GMRES is used for the main solver and the convergence threshold for GMRES is 10 −10 . The size of each subdomain is H = 1/2 and 64 processes are used. The second column (B Sp ) shows the CPU time for building an operator and the third column (S Sp ) shows the CPU time for one solve with S p . The fourth column (N I Sp ) shows the number of solves with S p that are required for convergence. The fifth column (TCPU) shows the total CPU time in seconds. The sixth to ninth columns show the corresponding results for S c . Note that the number of GMRES iterations required for convergence with S c is more than 3, which is not consistent with the spectral analysis done in [27]. Each solve time and building time is considerably higher for S c than S p , as expected. However, the number of GMRES iterations required for convergence is much higher for S p than S c . Additionally, the number of iterations for S p increases as φ decreases, while for S c it is bounded above. For this particular problem S p is a better choice for relatively large values of φ (i.e., φ = 10 −7 ), while S c is a better choice for any value of φ smaller than 10 −9 . 7. Conclusion. We have introduced a practical factorization of the Schur complement that arises from distributed optimal control of linear static systems. Due to the exact representation, if the range-space method is applicable, then one solve with the Schur complement is sufficient to obtain an optimal control solution. However, the Schur complement becomes ill-conditioned for large values of φ and as the mesh is refined. For example, the Schur complement that arises from a 3D cantilever problem with real material properties of rubber is prone to ill-conditioning when a relatively large regularization value is used. In such a case, an inaccurate solution is likely to be obtained if the range-space method is used, so solving a full KKT system simultaneously using a Krylov iterative method with a good preconditioner is recommended. The complex factorization of the Schur complement introduced in this paper can be used with a Schur-complement based preconditioner. The comparison between the approximate Schur complement of Pearson and Wathen and the complex factorization as a preconditioner in the full-space method shows promising results for the complex factorization even in the context of the full-space method. The Schur complement solve is done with the parallel domain decomposition linear solver FETI-DPH. The scalability of FETI-DPH (both numerical and parallel) as well as its dependence on φ are studied in two academic problems: a thermal 2D problem and a structural 3D problem. Edge-based rigid body modes augmentation is able to bring the number of iterations down, but further research is necessary to find an optimal augmentation in the context of PDE-constrained distributed optimal control. The combination of exact representation of the Schur complement and good scalability of the FETI-DPH solver in addition to extensibility of the representation indicates promise for use of the complex factorization in more complicated and practical problems.