Inner product free iterative solution and elimination methods for linear systems of a three-by-three block matrix form

Abstract

Large scale systems of algebraic equations are frequently solved by iterative solution methods, such as the conjugate gradient method for symmetric or a generalized conjugate gradient or generalized minimum residual method for nonsymmetric linear systems. In practice, to get an acceptable elapsed computing time when solving large scale problems, one shall use parallel computer platforms. However, such methods involve orthogonalization of search vectors which requires computation of many inner products and, hence, needs global communication of data, which will be costly in computer times. In this paper, we propose various inner product free methods, such as the Chebyshev acceleration method. We study the solution of linear systems arising from optimal control problems for PDEs, such as the edge element discretization of the time-periodic eddy current optimal control problem. Following a discretize-then-optimize scheme, the resulting linear system is of a three-by-three block matrix form. Various solution methods based on an approximate Schur complement and inner product free iterative solution methods for this linear system are analyzed and compared with an earlier used method for two-by-two block matrices with square blocks. The convergence properties and implementation details of the proposed methods are analyzed to show their effectiveness and practicality. Both serial and parallel numerical experiments are presented to further investigate the performance of the proposed methods compared with some other existing methods.

Introduction

Large scale linear systems, in particular for those arising from 3D partial differential equation problems considered here, are mostly solved by preconditioned Krylov subspace iterative methods. Such methods, of conjugate gradient (CG) type, minimum residual (MINRES [1]) or generalized minimum residual (GMRES [2]) type are based on computed search vectors, which must be orthogonalized and normalized, and involve hence inner products to compute the updated iterative solutions and search vectors used in the corresponding methods. When solving very large scale problems, it is preferable to use parallel computer facilities and split the computations in subdomains as in domain decomposition methods or subgraphs as in algebraic multilevel or similar methods, see e.g. [3], [4]. Then each processor is involved with local computations on subdomains or subgraphs. Krylov subspace method involves updates of solution vectors in the form $x≔x+\alpha y$, where $x$, $y$ are vectors and $\alpha$ is a scalar. Such computations can be done locally on subdomains and computer cores. However, the computation of the scalar coefficients and orthogonalization of the search vectors involve computation of global inner products, i.e., global information from every core and every node of the machine. This computation can be done on a single computer processor, but it needs communications of subvector parts and, when computed, global distributions of the inner products to all processors or cores. Alternatively, the global inner product can be computed on all processors, but the local vector information must then anyway be sent from each processor to all other processors, which clearly needs much communication of data. Such global communication of data can be very time consuming for solving problems of large sizes, and is a waste of time, since no, or only little, computation can be done during such communications. The lower bound of the latency of such an operation will grow with increasing parallel site because the parallelism is limited to tree-like structure with the number of leaves determined by the number of cores. Hence the scalar process will eventually become the limiting factor for the available parallelism in the algorithm. In future generations of parallel computations, global communication will cause even more costly global synchronization steps. There have been several efforts to hide such global communications. One can try to overlap them with other communications and calculations, which can include calculations done in the preconditioners. Practical algorithms can be used, see, e.g., [5], [6], but they involve then additional dot products and more matrix–vector products.

The classical Chebyshev iteration method is free of such inner products but to be effective requires accurate estimates of the lower and upper eigenvalue bounds. We present an efficient preconditioner where this is possible. As pointed out, for example, in [7], use of a Chebyshev iteration instead of conjugate gradient acceleration method can decrease the computer time significantly on a hypercube device. Therefore, it is of interest to consider inner product free iterative methods which do not require any global orthogonalization steps as an alternative to Krylov subspace methods. This is one of the topics of this paper.

We deal with algebraic linear systems arising, for instance, from optimal control problems for partial differential equations involving a state, a control and an adjoint, i.e., Lagrange multiplier variable; see, e.g., [8], [9], [10], [11], [12], [13] and the references therein. In particular, we deal with the eddy current electromagnetic optimal control problems with a time-periodic control function. Such systems lead to complex valued three-by-three block matrix forms. Due to the large size of the system, it is of practical importance to reduce it. We consider two elimination methods, one based on elimination of both the control and adjoint variables and one based on elimination of only one of them, which leads to a single equation but corresponding to a fourth order differential operator and a two-by-two block matrix with square matrix blocks, corresponding to second order differential operators, respectively.

For the first elimination approach, we use a preconditioner in the form of a product of two second order operators in discretized form. We also consider a parameterized variant of this preconditioner and the resulting iterative refinement methods. For the second elimination approach, we consider an approach based on first rewriting the equations, to some extent similar to an alternative direction type of block equation splitting (see, e.g., [14]). Therefore, the method is called alternating diagonal iteration sequential preconditioning methods, i.e. ADISP.

We mention that a somewhat related method is the Hermitian and skew-Hermitian splitting (HSS) method [15] for non-Hermitian positive definite linear systems, which is also an efficient inner product free method. This method is similar to the ADISP method. Some generalized forms of the HSS method have been proposed and analyzed for solving other type of problems. In particular, it has been modified to solve complex symmetric linear systems [16], [17] and distributed optimal control problems [18]. We show that it is not competitive with the other methods. During each iterative step in the ADISP method, we solve a linear system consisting of the sum of the two block matrices which appear in the two-by-two block matrix. One of the matrices is the complex conjugate of the other matrix. For this method, we derive a convergence factor estimate that holds uniformly in the finite element mesh parameter, the periodic angular frequency and the regularization parameter, i.e., the control cost parameter.

It should be noted that an earlier version of block matrix preconditioner has been used for time-harmonic parabolic and eddy current optimal control problems [8], [13]. This preconditioner has a practical implementation algorithm when used for Krylov subspace acceleration and it can lead to parameter independent approximation for the eigenvalue distribution of the preconditioned matrix. See also several earlier publications dealing with this special type of preconditioning [8], [13]. One of the earlier results dealt with equations of eddy current optimal control type, and those equations are considered here also. When accelerated with the GMRES method, the robust structured preconditioner in [12] outperforms the block-diagonal and block-triangular preconditioned MINRES and GMRES methods in [8], [19], [20]; see also the numerical results in [12], [21]. We remark that the preconditioner in [12], [21], [22] and those proposed in Section 3 can lead to accurate estimates of the upper and lower eigenvalue bounds and can hence efficiently be used for a Chebyshev iteration instead of the GMRES iteration method.

The so arising complex valued systems for the proposed methods can be solved in real arithmetics by rewriting the complex equation in two real valued forms. The then arising system has a two-by-two block structure such as we deal with here and can hence be solved by preconditioned Krylov subspace methods. Methods of this type have appeared in various publications, see e.g., [23], [24]. Hence, both methods lead to coupled inner–outer iteration methods and it is important to examine how the rate of convergence of the outer iteration process depends on the inner iteration stopping criterion.

The major aim of the paper is to present efficient outer iterative methods. However, as all methods deal with solving an inner system for a discretized partial differential equation, it is important to consider also efficient iterative methods for inner systems. For this purpose, some form of multigrid methods must be used. As pointed out in [25], [26], this is a more complicated issue since here we deal with Maxwell’s, i.e., curl-curl equation, than for the standard Poisson type equations. In this paper, we take advantages of the aggregation-based algebraic Multigrid (AMG) method [27], [28], [29], which however may not be the optimal choice. Presentation of more efficient inner solvers should henceforth also be further considered in the future.

The paper is composed as follows. In Section 2, we present the model problem on which the methods will be tested and the two elimination methods used. Section 3 contains a presentation of the preconditioner to be used for the fourth order discretized differential operator and the resulting iterative refinement method. Section 4 contains a derivation of the rate of convergence for the ADISP method which is then applied to this model problem. In Section 5, we consider the Chebyshev acceleration of the proposed iterative solution methods. Section 6 reviews some standard block preconditioners. Finally, numerical experiments are provided in Section 7.