1 Introduction
Many works have focused on efficient iterative methods for the solution of large linear systems with coefficient matrix in saddlepoint 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/Toeplitzlike structure, some saddlepoint 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 saddlepoint structure is that of constructing preconditioning for twobytwo 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 3step 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 matrixvalued symbol. In Section 2, we introduce an iterative 3step procedure that can be exploited for solving linear systems with structured coefficient matrices in saddlepoint 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 bilevel 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 saddlepoint 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
where
(1) 
are the Fourier coefficients of and
(2) 
is the Fourier series of . The function is called generating function of .
The type of domain (either onedimensional 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 matrixvalued 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 matrixvalued function .
Type of generating function  Associated Toeplitz matrix  

univariate scalar  unilevel scalar  
variate scalar  level scalar  
univariate matrixvalued  unilevel block  
variate matrixvalued  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 matrixvalued 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 of
generated by a Hermitian matrixvalued function.Theorem 1 ([serra1998, serra99])
Let be a Hermitian matrixvalued 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 blockToeplitz matrix generated by a matrixvalued 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 fullrank 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 presmoother and times a postsmoother 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 TwoGrid 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 Vcycle method and the choice defines the Wcycle 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.
Theorem 1
([RStub]) Let and be defined as in the TGM algorithm and assume that no presmoothing 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:
Then and
1.3 Multigrid method for block Toeplitz systems
If the coefficient matrix is of the form , , with being a matrixvalued 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
(3) 
or the relaxed Jacobi method with relaxation parameter
(4) 
where represents the Fourier coefficient of present on the block diagonal of , corresponding to the multiindex .
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 matrixvalued trigonometric polynomial, , such that there exist unique and such that
(5) 
The latter means that only one eigenvalue function of has exactly one zero in and is positive definite in . Define
as the eigenvector function of
associated with . Moreover, define . Let us define a matrixvalued trigonometric polynomial such that
(6) which implies that the trigonometric function
is welldefined for all .

(7) 
(8) where is a constant.
In order to have optimal convergence of the TGM it is sufficient compute a grid transfer operator as
(9) 
where is a multilevel cutting matrix obtained as the Kronecker product of proper unilevel cutting matrices and is a multilevel blockToeplitz matrix generated by .
2 A 3step procedure
Let be a nonsingular 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
(10) 
with symmetric positive definite (SPD), nonnegative 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 TwoGrid method following the “unknownbased” 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
(11) 
Then, we multiply both sides of (10) to the left by
(12) 
obtaining the system
(13) 
with
(14) 
where .
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 TwoGrid 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 “unknownbased” multigrid approach where we consider separately the diagonal blocks associated to the velocity and the pressure. Consequently, we can consider a block TwoGrid method with the prolongation operator given by
(15) 
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 postsmoothing. However, after one application of the grid transfer operator in (15), the matrix at the coarse level has the form
(16) 
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
(17) 
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 Vcycle 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
(18) 
with SPD, nonnegative 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
where
For the considerations that we made when we described our 3SP, for each we have that is SPD, is nonnegative definite, and .
3 Block multigrid methods for indefinite matrices
In this section we analyze and rewrite the threestep 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 threestep 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.
Proposition 1
Consider the iterated three steps procedure described in Section 2 applied to system (10). At level , define the prolongation and restriction operators as
(19) 
(20) 
After iterations of 3SP, the resulting coefficient matrix is given by
(21) 
Proof
The assertion follows immediately by direct computation.
4 Saddlepoint 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
(22) 
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
where
Higher order elements are needed for the discretization of the velocities than for the discretization of the pressure , e.g., Q2Q1 TaylorHoodelements. Alternatively, linear elements can be used for the displacement but using a subdivision of the element, corresponding to a 2times finer—or once refined—mesh for the displacements. We choose the latter approach that is often called Q1isoQ2/Q1 [MR1115205]. The properties of the system matrix are similar to those of the system matrix produced by Q2Q1 elements. The advantage of using Q1isoQ2/Q1 is that the discrete operator representing is more sparse than that for Q2Q1 and thus cheaper to apply. The resulting algebraic system of equations given by
(23) 
where and , so . In particular, the matrix is nonnegative and it is the permutation of a bilevel 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 bilevel Toeplitz matrixwhere
(24) 
with .
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 bilevel system. This property will be used later when we will exploit the information of the permuted matrix where is the bilevel block Toeplitz generated by . In particular,
(25) 
where
(26) 
4.2 Multigrid Method for the linear system
To solve efficiently the system (23) with the method proposed in Section 3 we investigate the structure of the coefficient matrix . It is a block matrix of the form
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 .
In order to apply the multigrid method of Section