Many works have focused on efficient iterative methods for the solution of large linear systems with coefficient matrix in saddle-point form [MR3564863, MR3439215, benzi2005numerical]. Indeed, many preconditioned Krylov subspace methods have been developed based on the study of the spectral properties of the coefficient matrix [MR3564863, MR2220678, benzi2006eigenvalues, locspsb1]. When we deal with multilevel and multilevel block systems the performance of some known preconditioners can deteriorate [NSV]. Multigrid methods are in general a good choice in terms of computational cost and optimal convergence rate. Indeed, they permit to achieve a fast convergence rate by constructing a proper sequence of linear systems of decreasing dimensions. Moreover, when the coefficient matrix has a particular multilevel (block) Toeplitz structure, the convergence analysis for such methods can be obtained in compact form exploiting the concept of symbols [FS1, FS2].
The knowledge and the use of the symbol has been successfully employed for designing very efficient preconditioning and multigrid techniques in various settings and several applications [MR3689933, DFFSS, MR4188673, sesana2013spectral]. However, even in the case where the coefficient matrix is somehow related to a Toeplitz/Toeplitz-like structure, some saddle-point problems are intrinsically indefinite. This limits the possible convergence results that can be exploited.
In many works [MR3564863, MR2220678, MR3689933] a common choice for problems with saddle-point structure is that of constructing preconditioning for two-by-two block systems involving the Schur complement. In many cases stemming from the approximation of PDE, even though the whole problem is indefinite, the latter procedure is based on the solution of a linear system with a positive definite coefficient matrix.
In this paper we consider the approach suggested by [MR3439215] that transforms the system in a new one that permits us to define an efficient 3-step procedure whose main body is formed by a symbol based block multigrid method.
In particular we construct the grid transfer operators and we focus on how to proceed for verifying the conditions which lead to convergence and optimality of the method.
The outline of the paper is as follows. In Subsection 1.1 we provide the notation and the necessary results on multigrid methods. We focus especially on the case where it is applied to Toepliz matrices with either scalar or matrix-valued symbol. In Section 2, we introduce an iterative 3-step procedure that can be exploited for solving linear systems with structured coefficient matrices in saddle-point form efficiently. Then, in Section 3 we show how the presented procedure can be studied as a block multigrid method and in Section 4 we explore the convergence results for the derived method when it is applied on the bi-level system stemming from the discretization of the Stokes equation in 2D. In particular, we focus on the structure and on the spectral features of the involved matrix sequences in order to prove the efficiency of the multigrid strategy with the proposed projection operators given in Section 3. Finally, we conclude with Section 5 presenting selected numerical tests that confirm the efficiency of the proposed multigrid method in comparison with known different strategies for saddle-point structured problems.
1.1 Notation and preliminary results
This Subsection is devoted to fix the notation and present the key concepts on multigrid methods and Toeplitz matrices.
We report the main definitions and properties of Toeplitz and circulant matrices, for simplicity, in the scalar unilevel setting. We only give the generalization in the block (or multilevel block) case of the results that will be exploited for the purpose of the paper.
An Toeplitz matrix is a matrix that has equal entries along each diagonal, and can be written as
Among the Toeplitz matrices, we focus on a particular case where the Toeplitz matrix is associated with a function belonging to and periodically extended to the whole real line. In this scalar and unilevel case the matrix is defined as
are the Fourier coefficients of and
is the Fourier series of . The function is called generating function of .
The type of domain (either one-dimensional or -dimensional ) and of the codomain (either the complex field or the linear space of complex matrices) of gives rise to different kinds of Toeplitz matrices, see Table 1 for an overview. We will denote the generating function by bold when it is a matrix-valued function. When we want to highlight that the generating function is multivariate, we specify its argument in bold, where . The structure and the size of the Toeplitz matrix is clearly related to its generating function. In the most general block multilevel case, given an -index we indicate by the size of the -level block Toeplitz matrix, denoted by , generated by an -variate matrix-valued function .
|Type of generating function||Associated Toeplitz matrix|
|univariate scalar||unilevel scalar|
|-variate scalar||-level scalar|
|univariate matrix-valued||unilevel block|
|-variate matrix-valued||-level block|
In the whole paper we adopt the following notation on norms and relations between positive definite matrices. If is a positive definite matrix, denotes the Euclidean norm weighted by on . If and are Hermitian matrices, then the notation means that is a nonnegative definite matrix. Finally, given a matrix-valued function defined on the cube we denote by .
The generating function plays an important role for investigating some analytic properties of the related Toeplitz matrix. Under particular hypotheses, can provide important spectral information on . In this case we say that is the spectral symbol of the matrix sequence
. In the following, we report an important localization result on the eigenvalues ofgenerated by a Hermitian matrix-valued function.
Theorem 1 ([serra1998, serra99])
Let be a Hermitian matrix-valued function with . Let be the essential infimum of the minimal eigenvalue of , be the essential supremum of the minimal eigenvalue of , be the essential infimum of the maximal eigenvalue of , and be the essential supremum of the maximal eigenvalue of .
If then all the eigenvalues of belong to the open set for every . If then all the eigenvalues of belong to the open set for every .
If and is the unique zero of such that there exist positive constants for which
then the minimal eigenvalue of goes to zero as .
1.2 Multigrid methods
In this subsection we briefly report the relevant results concerning the convergence theory of algebraic multigrid methods [RStub] especially when the coefficient matrix is a block-Toeplitz matrix generated by a matrix-valued trigonometric polynomial [block_multigrid2].
In the general case, we are interested in solving a linear system where is a Hermitian positive definite matrix. Assume and define two sequences of full-rank rectangular matrices and . Let us consider the stationary iterative methods , with iteration matrix , and , with iteration matrix .
An iteration of an algebraic Multigrid Method (MGM) is given by the following algorithm.
The steps in lines and consist, respectively, in applying times a pre-smoother and times a post-smoother of the given iterative methods. The steps in lines to define the “coarse grid correction”. When and are equal to 1 the Algorithm 1 represents one iteration of the Two-Grid method (TGM). In this case a linear system with coefficient matrix has to be solved. Consequently, if the coefficient matrix has a large size, should be taken such that becomes small enough for solving cheaply step 1. In the latter situation, the choice defines the V-cycle method and the choice defines the W-cycle method.
The next Theorem provides the conditions which are sufficient for the convergence of the TGM in the particular case where the matrices at the coarse levels are computed according to the Galerkin approach (). Moreover, the theorem guarantees the linear convergence of the method, which means that the number of iterations needed by the algorithm to reach a given accuracy is bounded from above by a constant independent of the size of the problem. In this case we will say that the method is convergent and optimal.
([RStub]) Let and be defined as in the TGM algorithm and assume that no pre-smoothing is applied . Let us define
the iterative matrix of the TGM. Assume that
the smoother fulfills the smoothing property:
The grid transfer operator fulfills the approximation property:
1.3 Multigrid method for block Toeplitz systems
If the coefficient matrix is of the form , , with being a matrix-valued trigonometric polynomial, , sufficient conditions for the convergence and optimality of the TGM can be expressed in relation to . In particular, in [DFFSS] the authors show that validating the smoothing condition (a) of Theorem 1 for is possible, for example, when considering as smoother either the relaxed Richardson method with relaxation parameter
or the relaxed Jacobi method with relaxation parameter
where represents the Fourier coefficient of present on the block diagonal of , corresponding to the multi-index .
In order to validate the approximation property (b) of Theorem 1, in [block_multigrid2] sufficient conditions that should be fulfilled by the grid transfer operators are derived. In the following we report a brief description of the result [block_multigrid2] in the general multilevel block setting.
Suppose is a matrix-valued trigonometric polynomial, , such that there exist unique and such that
The latter means that only one eigenvalue function of has exactly one zero in and is positive definite in . Define
as the eigenvector function ofassociated with . Moreover, define . Let us define a matrix-valued trigonometric polynomial such that
which implies that the trigonometric function
is well-defined for all .
where is a constant.
In order to have optimal convergence of the TGM it is sufficient compute a grid transfer operator as
where is a multilevel cutting matrix obtained as the Kronecker product of proper unilevel cutting matrices and is a multilevel block-Toeplitz matrix generated by .
2 A 3-step procedure
Let be a non-singular matrix. Assume that we want to solve a sequence of linear systems of increasing dimension with as coefficient matrix, with size of the form . Moreover, assume that, for a fixed , we have a system of the form
with symmetric positive definite (SPD), non-negative definite, and . To solve efficiently the system (10) we define an iterative procedure consisting of 3 steps (3SP):
transform the system into a new one by making use of two invertible matrices and .
Apply one step of standard Two-Grid method following the “unknown-based” multigrid approach.
“Symmetrize” the problem with a proper variable transformation.
Concerning the first step, in particular, we transform the system following the approach suggested by [MR3439215].
In particular, we fix a positive parameter , e.g., where . We perform the change of variables
Then, we multiply both sides of (10) to the left by
obtaining the system
Note that once we solve (13), the solution of the system (10) is given by relation (11). Moreover, we highlight that the computational cost with the proposed approach is moderate compared to that of the original system thanks to the implementation described in [MR3439215].
After step one, the transformed system is not symmetric. However, the block is a SPD matrix as long as the parameter is such that . Equivalently, the matrix in block position is SPD, since it is given by .
The next step consists in applying a Two-Grid method to solve system (13) where the coefficient matrix at the finer level is given by the matrix in (14). For this purpose we consider one step of the “unknown-based” multigrid approach where we consider separately the diagonal blocks associated to the velocity and the pressure. Consequently, we can consider a block Two-Grid method with the prolongation operator given by
where and are the prolongation operators chosen for solving efficiently the scalar systems with coefficient matrix and , respectively. Moreover, the method should be completed by a basic scheme using only a single step of damped Jacobi post-smoothing. However, after one application of the grid transfer operator in (15), the matrix at the coarse level has the form
which, apart from the minus in block position , is a matrix of a form which resembles that of the system (10). Then, we apply the third step of our procedure with consists in the variable transformation obtained multiplying both (16) and the right hand side on the left by the matrix
This last step restores the symmetry of the matrix and gives to (16) the same block form of the original coefficient matrix in (10). Then, we can apply recursively our three step procedure with updated and and the proper in order to define a V-cycle method.
2.1 Recursive definition of matrices
If we have a matrix of the form in equation 10, our three step procedure suggests to define a sequence of matrices with a recursive definition.
The base of the recursive definition is given by the formula
with SPD, non-negative definite, and .
At level , suppose that and are the prolongation operators chosen for solving efficiently the scalar systems with coefficient matrix and , respectively, where is such that .
Then, the inductive step of our recursive definition is given by
For the considerations that we made when we described our 3SP, for each we have that is SPD, is non-negative definite, and .
3 Block multigrid methods for indefinite matrices
In this section we analyze and re-write the three-step procedure of the previous section in order to develop one block multigrid method for the problem in (10). Precisely we answer the following main questions.
Can the three-step procedure be seen globally as one main block multigrid procedure for the problem in (10)?
If yes, which are the sequences of grid transfer operators and , needed for the construction of the coarse coefficient matrices?
The following proposition, answers to 1. and 2., analyzing how the structure of the matrices changes after applying times the 3SP.
After iterations of 3SP, the resulting coefficient matrix is given by
The assertion follows immediately by direct computation.
4 Saddle-point matrices stemming from the Stokes equation
In this section we are interested in applying the multigrid method described before to large linear systems stemming from the finite element approximation of the Stokes equation. The Stokes equation has the form
where is the velocity and is the pressure, by we denote the Dirichlet boundary part of and by the Neumann part. For the construction and analysis of the system matrices we consider the domain , homogenous Dirichlet boundary conditions are attained for and , homogenous Neumann boundary conditions for and .
4.1 Finite element discretization
The weak form of (22) is given by
Higher order elements are needed for the discretization of the velocities than for the discretization of the pressure , e.g., Q2-Q1 Taylor-Hood-elements. Alternatively, linear elements can be used for the displacement but using a subdivision of the element, corresponding to a 2-times finer—or once refined—mesh for the displacements. We choose the latter approach that is often called Q1-iso-Q2/Q1 [MR1115205]. The properties of the system matrix are similar to those of the system matrix produced by Q2-Q1 elements. The advantage of using Q1-iso-Q2/Q1 is that the discrete operator representing is more sparse than that for Q2-Q1 and thus cheaper to apply. The resulting algebraic system of equations given by
where and , so . In particular, the matrix is non-negative and it is the permutation of a bi-level Toeplitz matrix with partial dimension . Precisely, let be the
th column of the identity matrix of size, we can define a proper permutation matrix, , , such that the th column of , is . The matrix transforms the matrix into the following a bi-level Toeplitz matrix
We consider an analogous permutation matrix, , , such that the th column of , is , with being the th column of the identity matrix of size . We can transform the global matrix in a bi-level system. This property will be used later when we will exploit the information of the permuted matrix where is the bi-level block Toeplitz generated by . In particular,
4.2 Multigrid Method for the linear system
with , and . Hence, it can be seen as the matrix at the level 0 in the formula (18), with . Because of the structure of , . Then, when we choose the positive parameter , we can write .