Refined saddle-point preconditioners for discretized Stokes problems

This paper is concerned with the implementation of efficient solution algorithms for elliptic problems with constraints. We establish theory which shows that including a simple scaling within well-established block diagonal preconditioners for Stokes problems can result in significantly faster convergence when applying the preconditioned MINRES method. The codes used in the numerical studies are available online.

together with suitable (Dirichlet, Neumann or mixed) boundary conditions. Stokes problems typically arise when modelling the flow of a slow-moving fluid such as magma in the Earth's mantle, see [18]. In our setting v denotes the flow velocity, p is the pressure, and f represents a source term that drives the PDE system. The associated boundary value problem is usually posed on a bounded domain ⊂ Rd , d ∈ {2, 3}. Stokes problems also arise in a natural way when the (unsteady) Navier-Stokes equations are simplified using classical operator splitting techniques, see [6].
We suppose that the boundary value problem is discretized using standard mixed finite elements. That is we take {φ i } i=1,...,n v as the finite element basis functions for the velocity components (we assume that the same approximation space is used for each one), and {ψ i } i=1,...,m for the pressure; so that n v and m are the number of velocity and pressure grid nodes respectively. Having set up the associated velocity basis set (for example, { φ 1 , . . . , φ 2n v } := {(φ 1 , 0) T , . . . , (φ n v , 0) T , (0, φ 1 ) T , . . . , (0, φ n v ) T } in two dimensions), the resulting discrete Stokes system is the saddle-point system, where A ∈ R n×n (with n =dn v ) is the vector-Laplacian matrix given by and B ∈ R m×n is the divergence matrix The vectors v, p are discretized representations of v, p, with f , g taking into account the source term f as well as nonhomogeneous boundary conditions. The matrix C is the zero matrix when a stable finite element discretization (such as the Q 2 -Q 1 Taylor-Hood element) is used, and is the stabilization matrix otherwise. We assume that A is symmetric positive definite, which is the case when a Dirichlet condition is imposed on at least part of the boundary. The matrix C is always positive semidefinite. For consistency with the continuous Stokes system the matrix B should satisfy 1 ∈ null(B T ) in the case of enclosed flow (see, e.g., [8,Chapter 3]). However, other vectors may also lie in the nullspace of B; these are artefacts of the discretization, or arise from the imposition of essential boundary conditions. The matrix system (1.1) is of classical saddle-point form. 1 There has been a great deal of research devoted to solving systems of the form (1.1) using preconditioned iterative methods; see [2] for a definitive review. This body of work is relevant to any linear system that is generated by a mixed approximation; see [4,Chapter 3] for a characterization. To state the key spectral properties, it is useful to let where A ∈ R n×n is symmetric positive definite as above, C ∈ R m×m is symmetric positive semidefinite, B ∈ R m×n with m ≤ n and rank(B) = r ≤ m. We suppose that the (negative) Schur complement of A, has rank p. Then under these conditions A has n positive eigenvalues, p negative eigenvalues and m − p zero eigenvalues [2, page 21]. A widely studied block diagonal preconditioner for A is given by where H ∈ R n×n is some symmetric positive definite approximation to the Schur complement S. In the case where H = S and C = 0 (whereby B must be full rank for S, and hence P 1 , to be invertible), it is known that the eigenvalues of the preconditioned system are given by [14,17] λ(P −1 1 A) ∈ 1, (1.4) and in the case where the approximation of S (or indeed A) is inexact the preconditioner is frequently found to be extremely effective also. When the condition on C is weakened to allow the matrix to be symmetric positive semi-definite, it can be shown that 2 λ(P −1 1 A) ∈ −1, . (1.5) We note that, for the results (1.4) and (1.5), we have assumed invertibility of S in order for P 1 itself to be invertible, as P 1 includes the exact Schur complement. In the remainder of this paper, however, we consider situations where the Schur complement could be singular, but construct an inexact approximation H which is invertible.
In the specific case of the Stokes equations, the approximate Schur complement is either the mass matrix associated with the pressure approximation space 3 or an approximation. Common approximations of Q are its diagonal (see [23,27]), a lumped version (see [24]), or a Chebyshev semi-iteration method applied to Q (see [12,13,30]). We will study a refined version of the classical preconditioner in this work: instead of taking S ≈ H , our idea is to incorporate a scaling constant α > 0 and investigate using as a potential preconditioner for A. Intuitively there is little reason to assume that the matrix P α would be a more effective preconditioner than P 1 : by scaling the Schur complement we are after all moving the preconditioner 'further' from the ideal preconditioner P 1 . Remarkably, however, we frequently observe a significant improvement in the Stokes case. This improvement is justified theoretically herein. We also explain why setting a large value of α can significantly improve the performance of the iterative solver when a stabilized mixed approximation is employed.
We highlight a related discussion on [23, page 1361] where a small scaling parameter was considered: the motivation for this was to reflect the behavior of the stabilization parameter β multiplying C within the Schur complement approximation. May and Moresi [16] scaled H = Q by the (fixed) viscosity of the fluid; the same scaling is applied in the Cahouet and Chabard preconditioner for generalized Stokes problems [5]. However, none of these investigate the optimal choice of scaling parameter.
Before we continue, let us fix notation. We order the eigenvalues of P −1 α A from smallest to largest, so that where p = rank(S) ≤ m and S is in (1.2). Additionally, the notation (F, G) is used to denote the generalized eigenvalue problem Fv = λGv. The n × n matrix formed by extracting the diagonal of F ∈ R n×n will be denoted diag(F).

Spectral equivalence bounds
Extensions to existing eigenvalue bounds for the Stokes problem are discussed in this section. We analyse the "ideal" Stokes preconditioner (1.6) first, but we also discuss bounds for efficient "inexact" variants. These results provide informal motivation for modifying the standard saddle-point preconditioner for the Stokes equations. Refined eigenvalue estimates applicable in a Stokes setting are presented in Sect. 3.

General saddle-point systems
We first wish to fix ideas using a general saddle-point system A, preconditioned by P α . We characterize the eigenvalues of P −1 α A using the following theorem. Although the result is simple, and is similar in flavour to results in many other papers (e.g., [3,11,19,23]), it forms the basis of our analysis and so we provide a proof for completeness.
We highlight that this corresponds to an exact application of the (1, 1)-block A within the preconditioner.

Theorem 2.1 Consider the generalized eigenvalue problem
with A ∈ R n×n symmetric positive definite, C ∈ R m×m symmetric positive semidefinite and B ∈ R m×n , m < n. Assume that rank(B) = r ≤ m. Then, where λ ∈ R, x ∈ R n and y ∈ R m , with x and y not simultaneously zero vectors. If C = 0, then Case II occurs if and only if λ = 0.
Proof Equation (2.1) is equivalent to We consider Cases I-III separately.
Case I If λ = 1 then (2.2) implies that B T y = 0, so either y = 0 or y ∈ null(B T ), y = 0. If y = 0 then (2.3) implies that B x = 0, so that x ∈ null(B). There are n − r linearly independent such vectors. Otherwise, y ∈ null(B T ) with y = 0. However, premultiplying (2.3) by y T then gives that α y T H y = −y T C y. Since H is positive definite, C is semidefinite and α > 0, this cannot hold. Thus, if λ = 1 then y = 0. On the other hand, if y = 0 we know from (2.2) that λ = 1, since x = 0 and A is positive definite, so λ = 1 if and only if y = 0. Accordingly, 1 is an eigenvalue of (A, P α ) with multiplicity n − r and eigenvectors [x T , 0 T ] T , x ∈ null(B).
Case II We now assume that y ∈ null(B T ), y = 0. From Case I we know this implies that λ = 1. Then, (2.2) shows that x = 0. From (2.3) it follows that λ and y satisfy the generalized eigenvalue problem −C y = λα H y. Thus y must simultaneously be an eigenvector of (−C, αH ) and in the nullspace of B T . At most m − r linearly independent vectors satisfy this requirement. If C = 0, this case only arises if λ = 0.
On the other hand, if λ = 0 and C = 0 then (2.2) and (2.3) imply that x = 0 and y ∈ null(B T ), so that Case II applies.
Case III Otherwise, we know that λ = 1, x = 0, y / ∈ null(B T ). We can rearrange (2.2) for x and substitute into (2.3) to give We see that it is possible to describe the eigenvalues of P −1 α A in terms of A, B, C, H and α. We also note that when C = 0 (as arises when solving the Stokes equations using stable finite elements), Case III describes all eigenvalues not equal to 0 or 1. This is a good place to pause to consider the implications of Theorem 2.1 and the effect of scaling P α on the eigenvalues of the preconditioned matrix for the Stokes equations. Trivially, eigenvalues satisfying Case I are positive (since λ = 1) while any eigenvalues satisfying Case II are non-positive, since C is semidefinite and H is positive definite. The remaining eigenvalues of P −1 α A may be positive, negative or zero and the inertia of P −1 α A must be the same as that of A. However, because C is semidefinite and A and H are positive definite, any positive eigenvalue approaches 1 from above as α increases. On the other hand, we see that negative eigenvalues may approach zero from below as α increases, which can have a detrimental affect on the speed of convergence of preconditioned MINRES. For this reason, it is interesting and important to examine in greater detail the effect of α on the eigenvalues of P −1 α A.

The Stokes equations
The contributions we provide in this paper rely on the particular numerical properties of the Stokes equations, so we now present a framework for this problem by considering suitable approximations of the (1, 1)-block and associated Schur complement. We first note that, in practice, an effective preconditioner will not invert the (1, 1)block exactly as this will be very expensive computationally. However it is reasonable to assume, as in [23], that an approximation A may be constructed such that for some function g of the mesh parameter h. Applying a tailored multigrid method to approximate the action of A −1 , for example, will achieve this property with g(h) bounded away from zero independently of h. For stable finite element discretizations of the Stokes equation there exists an inf-sup constant γ , and a constant resulting from the boundedness of B, such that For an unstable discretization only the upper bound holds, and a lower bound is assumed as follows (provided p / ∈ span{1} in the case of enclosed flow): Furthermore we assume that there exist mesh-independent constants θ , guaranteeing the spectral equivalence of Q and the Schur complement approximation H , that is: Finally we use the boundedness of C to write for some mesh-independent constant . The properties assumed above all hold for the discretizations and approximations we use in this work. We are now in a position to recall Theorem 2.2 of [23] which in turn provides a bound for the convergence of preconditioned MINRES, see [29, The asymptotic convergence rate of preconditioned MINRES, denoted by e k , satisfies (2.10) We refer to [29] for discussion of the asymptotic convergence rate for such problems, and to [26,Chapter 3.2] for a definition and motivation of this quantity. We highlight at this stage that g(h) is solely determined by how one defines the (1,1) block in the preconditioner, for example using a multigrid method. As methods for approximating such matrices are very well-established, we will therefore regard this as a fixed number. From (2.10) therefore, we observe that the quantity controlling the 'average' convergence of the method is That is, if R 1 is maximized, then the 'best' average convergence is achieved. We note that = 0 when no stabilization is applied.
Remark In [23, Corollary 2.1], the authors proceed to demonstrate that if A is spectrally equivalent to A (i.e. with no dependence on h), then the convergence rate of the iterative scheme is independent of the mesh. The result (2.10), which assumes some dependence on h through the assumption (2.5), does however remain a highly valuable tool to analyse the consequences of parameter changes on the convergence rate, so it is important to state it here.
We should also highlight the fact that the function g(h) will not depend on h in practical applications (since spectrally robust methods such as multigrid can be used to approximate A). Therefore the lower bound in assumption (2.5) would ideally be replaced by a mesh-independent constant, C l say. To make this assumption useful in our setting, we may take g(h) = min{C l h δ , C l } for instance, where δ > 0 is sufficiently small that h δ remains roughly 1 for all h tested, whereupon (2.5) is satisfied and the dependence on h in (2.10) is nullified for all practical computations.

The effect of scaling
We now consider the result of applying the scaled preconditioner given by on this known theoretical result. Then within our assumptions (2.8), (2.9) for Theorem 2.2, we must replace θ 2 , 2 , with θ 2 /α, 2 /α, /α, in which case the asymptotic convergence rate becomes We now examine the behavior of R α as α ↑ ∞, starting with the case where a stable discretization is used (i.e. = 0). In this case In the case where a stabilized mixed method is used (i.e. = 0), we have The above discussion indicates that, for both stable and stabilized discretizations, it may be highly advantageous to increase the scaling parameter α in P α . In particular, increasing α nullifies the effect of the parameter in the expression for the average convergence rate. As α ↑ ∞, the predicted rate tends to 1 − g(h) 3/4 √ γ θ/ . Of course this argument is a heuristic, as we do not know from this working how large α must be to result in substantially faster convergence. While scaling a saddle-point preconditioner is a strategy that is commonly adopted by practitioners to accelerate convergence, tuning is invariably done without theoretical justification. We would like to fix this in the Stokes flow context. We provide justification for setting α to a moderately large value in Sects. 3 and 4, and the performance gains that are achievable when choosing a sensible scaling parameter are demonstrated in Sect. 5.

Refined estimates for the negative eigenvalues of Stokes problems
The bounds in the previous section suggest that large values of α in P α will reduce the condition number of P −1 α A and hence improve the convergence rate of preconditioned MINRES applied to Stokes problems. Fast convergence of Krylov subspace methods for symmetric indefinite problems is often attributed to nicely distributed eigenvalues, with clustered eigenvalues often sought. 4 Recalling the remarks after Theorem 2.1, we find that positive eigenvalues of P −1 α A cluster near one as α increases. Negative eigenvalues also cluster as α increases, but move towards the origin, which can delay the convergence of Krylov subspace methods.
Accordingly, it is instructive to more precisely characterize λ p , the negative eigenvalue of P −1 α A nearest the origin, for Stokes problems discretized by different finite elements. In particular, we examine stable Q 2 -Q 1 elements, and the two stabilized elements available in the Incompressible Flow & Iterative Solver Software (IFISS) [9,10,22] software. These are Q 1 -Q 1 elements with the stabilization approach of Dohrmann and Bochev [7] (see also [8,Sect. 3.3.2]) and Q 1 -P 0 elements stabilized as in [8,Section 3.3.2]. Note that for these Q 1 -P 0 elements the pressure mass matrix is diagonal, so that Q = diag(Q). We assume in this section that the (1, 1) block of P α is A and the (2, 2) block is either the pressure mass matrix or its diagonal.
To motivate our analysis, we compute the extreme negative and positive eigenvalues λ 1 , λ p , λ m+1 and λ m+n of P −1 α A as α varies, for a cavity problem discretized by Q 1 -Q 1 , Q 1 -P 0 and Q 2 -Q 1 elements. This is a widely considered problem in Stokes flow, which we define on = [−1, 1] 2 , with f = 0 and boundary conditions given by Since the flow is enclosed, the preconditioned system is singular with a single zero eigenvalue that is associated with a zero velocity and a constant pressure vector.
Tables 1 and 2 are consistent with the asymptotic results for large α in Sect. 2. We also note that λ p approaches the origin algebraically as α is increased. Other interesting trends also emerge. One intriguing feature of Q 1 -Q 1 elements is that, when H in P α is the diagonal of the pressure mass matrix, the eigenvalue λ p seems to be −0.25/α. On the other hand, when H is the full pressure mass matrix and α is large the eigenvalue λ p is almost (although not exactly) the same for all three element types.
Our next task is to develop good bounds for λ p and explain some of the phenomena we observe, so that we might choose a value of α that results in fast convergence of Krylov methods applied to Stokes problems. To do this we examine both Case II and Case III eigenvalues from Theorem 2.1.

Case III eigenvalues
We start by studying Lemma 2.3 in [23] (adapted to include α), which can also be obtained by bounding the Case III eigenvalues in Theorem 2.1. Importantly, when H = Q, the pressure mass matrix, the bound is remarkably tight. Although the same bound can be applied when H = diag(Q), we will see that the results are less informative since θ 2 γ 2 is far from ν min , the smallest value of ν in Theorem 2.1.   Table 1 Computed extreme eigenvalues of P −1 α A for the cavity problem, a mesh parameter of 2 −5 and H = diag(Q), the diagonal of the pressure mass matrix Since Q 1 -Q 1 and Q 1 -P 0 elements satisfy the ideal stabilization property (see [8,Sect. 3.3.2]), the largest eigenvalue of Q −1 C is less than or equal to 1 for these elements. Additionally, for stable Q 2 -Q 1 elements C = 0 so that (2.9) is trivially satisfied. Thus, for all three elements we can apply Lemma 3.1. Moreover, γ is bounded away from zero by a constant that depends on the element type but not on the mesh parameter h [8,Sect. 3.5]. The smallest value, ν min , of ν in Theorem 2.1 approximates γ 2 and is given in Table 3 for the cavity problem; for the obstacle problem this is introduced in Sect. 5. Tables 4 and 5 show the bound (3.1) and corresponding value of λ p for the cavity problem. We see that the bound is pessimistic when H = diag(Q) (with the exception of Q 1 -P 0 elements for which diag(Q) = Q). However, the bound is very accurate for all three elements when the full pressure mass matrix is used in P α . Additionally, when H = Q, it appears that the eigenvalue λ p is determined mainly by ν min , which varies only mildly between the different element types, and which is bounded away from zero independently of h. Qualitatively similar results are observed for the obstacle problem described in Sect. 5. Importantly, it seems that when H is the pressure mass matrix we can accurately bound λ p as α increases, which allows us to control the magnitude of this eigenvalue.   Table 4 Largest negative eigenvalue (λ p ) of P −1 α A, and bound (3.1) for the cavity problem, a mesh parameter of 2 −5 and H = diag(Q), the diagonal of the pressure mass matrix

Case II eigenvalues
Although the eigenvalue bound (3.1) is descriptive when we use the full pressure mass matrix in P α , it is rather pessimistic when we use its diagonal instead (except for Q 1 -P 0 elements). It would be useful to have an alternative means of quantifying λ p , when H = diag(Q), for Q 1 -Q 1 and Q 2 -Q 1 elements. The latter case appears to be difficult, since there are no Case II eigenvalues in Theorem 2.1. However, we see from Table 1 that for Q 1 -Q 1 elements λ p behaves like −0.25/α. We show in the rest of this section that this is indeed the case, and that this eigenvalue is associated with Case II in Theorem 2.1. Since it is possible to characterize the Case II eigenvalues for the full pressure mass matrix, and for Q 1 -P 0 elements, we extend our analysis to these cases for completeness. Case II eigenvalues satisfy −C y = λα H y, y ∈ null(B T ). Our approach for this analysis is to propose a basis for null(B T ), and then determine whether these basis Table 5 Largest negative eigenvalue (λ (1, n y + 1) (n x + 1, n y + 1) 1 2 3 4 Fig. 1 Diagram of mesh and nodes (left), and node numbering within each element (right) vectors are eigenvectors of the generalized problem (−C, αH ). To do so we require certain notation, and details of the finite element assembly process, that we describe here. We assume that there are n x elements in the x direction and n y elements in the y direction, so that the total number of elements is n el = n x n y (see Fig. 1). Although we restrict our attention to rectangular domains for simplicity, the same methodology can be used to analyse more complicated domains, as we discuss at the end of this section.
Let us first briefly introduce some notation to describe the assembly process. Let C k ∈ R × , Q k ∈ R × and diag(Q k ) ∈ R × , k = 1, . . . , n el , be the element matrices that are assembled to form C, Q and diag(Q). Additionally, let L ∈ R N ×m be the connectivity matrix that maps local pressure degrees of freedom on element k to the global pressure degrees of freedom 1, . . . , m, where N = n el . Then (3.2)

Q 1 -Q 1 elements
We now examine Case II eigenvalues for Q 1 -Q 1 elements. Let us begin by specifying the Q 1 -Q 1 connectivity matrix.

Lemma 3.2
Let the Stokes equations be discretized by Q 1 -Q 1 elements on a rectangular domain of square elements, with n x elements in one direction and n y elements in the other. Let L ∈ R N ×m be the Q 1 -Q 1 connectivity matrix that maps the N = 4n x n y local pressure degrees of freedom to the m global pressure degrees of freedom. Con-sider the (i, j)th node in the finite element mesh, where the node number is as in Fig. 1.
Then the corresponding column of L is given by 4 +e j ⊗ e i−1 ⊗ e 2 + e i ⊗ e 1 i = 2, . . . , n x , j = 2, . . . , n y , e j−1 ⊗ e n x ⊗ e 3 + e j ⊗ e n x ⊗ e 2 i = n x + 1, j = 2, . . . , n y , e n y ⊗ e 1 ⊗ e 4 i = 1, j = n y + 1, e n y ⊗ e i−1 ⊗ e 3 + e i ⊗ e 4 i = 2, . . . , n k , j = n y + 1, e n y ⊗ e n x ⊗ e 3 i = n x + 1, j = n y + 1, . . , n x + 1 and j = 1, . . . , n y + 1. In each Kronecker product e j ⊗ e i ⊗ e s , the vectors e j ∈ R n y , e i ∈ R n x and e s ∈ R 4 are the jth, ith and sth unit vectors of the appropriate dimension.
Since L has one element per row, that is, the connectivity matrix maps the constant vector to one of larger dimension (cf. Lemma 3.4 below).
As stated at the start of this section, we employ the stabilization approach of Dohrmann and Bochev, who define the stabilization matrix on the kth element to be where q = [ 1 4 , 1 4 , 1 4 , 1 4 ] T and | k | is the element area. With this choice it follows from (3.2) that the stabilization matrix C satisfies C = Q − | k |ww T , where w = 1 4 L T 1 N . It is straightforward to compute that null(C k ) = span{1 4 }. Since the connectivity matrix maps the constant vector to one of larger dimension (see (3.4)), null(C) = span{1 m }.
To show that we have a Case II eigenvalue, we must be able to show that ±1 is an eigenvector of −C y = λα diag(Q) y or, equivalently, of

Lemma 3.3 Let the Stokes equations be discretized by Q 1 -Q 1 elements on a rectangular domain of square elements. The eigenpairs of
Proof The result is obtained by straightforward computation.

, to vectors of length m via
where with k ∈ {−1, 1}, k = 1, . . . , n el . To proceed we require a technical result that shows that the v s lie in range(L). Proof We must be able to combine columns p of L in (3.3) to get v s , s = 1, . . . , 4. It is straightforward, although rather cumbersome, to show that which proves the result.

Since L(L T L) −1 L T is an orthogonal projector onto range(L), a consequence of Lemma 3.4 is that
(3.9) Importantly, this means that Lv 2 = v 2 = ±1 ∈ null(B T ). The final step is to combine (3.9) with Lemma 3.3 to show that v 2 is indeed an eigenvector of (−C, α diag(Q)), with corresponding eigenvalue −0.25/α. Proof From (3.2) we have that Cv s + λ s α diag(Q)v s = L T (C k + λ s α diag(Q k ))Lv s . Thus, using (3.8), (3.9) and Lemma 3.3, we find that which shows that (λ s , v s ) are eigenpairs of (−C, α diag(Q)).
Now we are in a position to determine Case II eigenvalues.

Theorem 3.6 Let the Stokes equations in two dimensions be discretized on a rectangular domain with square elements by Q 1 -Q 1 elements, and let H = diag(Q).
Then −0.25/α is the largest negative Case II eigenvalue of P −1 α A. Proof The vectors v s , s = 1, . . . , 4 are candidates for y in Case II eigenvalues in Theorem 2.1. As discussed above v 2 lies in null(B T ) before B is modified to accommodate any essential boundary conditions [20,21]. Thus −0.25/α is an eigenvalue of P −1 α A.
Our final task is to show that −0.25/α is the largest negative Case II eigenvalue. To do so we employ the approach of Wathen [28]. We let Q k represent the diagonal matrix whose diagonal entries are those of Q k to simplify notation. Since 1 m = L1 N , N = 4n el , any nonzero Case II eigenvalue λ must satisfy The eigenvalue λ = −0.25/α is thus precisely λ p that we observe in Table 1. Of course, certain boundary conditions may increase the dimension of null(B T ), in which case v 1 , v 3 and/or v 4 may lie in the nullspace of this modified matrix. For example, for the channel problem all four vectors lie in the nullspace.
For completeness we now consider the case where H = Q, the pressure mass matrix, which can be analysed in a very similar manner to H = diag(Q) above. As in the proof of Theorem 3.6, before B is modified to accommodate boundary conditions, null(B T ) = span{v 2 } and the eigenvalue −1/α is guaranteed to be a Case II eigenvalue of P −1 α A. A similar generalized Rayleigh-quotient analysis to that in the proof of Lemma 3.6 shows that this is the most negative Case II eigenvalue.
Again, v 1 , v 3 and/or v 4 may lie in the nullspace of B after essential boundary conditions are imposed. We stress that the difference between this and the diagonal pressure mass matrix approximation is that λ p = −1/α. Instead λ p is as in Table 2, and is bounded by (3.1).
More generally, we can bound the largest Case II eigenvalues for any preconditioner that is spectrally equivalent to Q, i.e., any preconditioner for which (2.8) holds. (2.8). Then the largest nonzero Case II eigenvalue of P −1 α A is bounded above by −θ 2 /α.

Corollary 3.8 Let the Stokes equations in two dimensions be discretized on a rectangular domain with square elements by Q 1 -Q 1 elements, and let H satisfy
Proof Any nonzero Case II eigenvalue satisfies where we have used Theorem 3.7.
Corollary 3.8 allows us to bound Case II eigenvalues for general preconditioners. Moreover, it gives some insight into whether it is likely that λ p is a Case II eigenvalue or a Case III eigenvalue.

Q 1 -P 0 elements
We now turn our attention to Q 1 -P 0 elements, which have one pressure degree of freedom per element, located at the element centre. A consequence is that the pressure mass matrix is diagonal, so that Q = diag(Q) = | k |I , where | k | is the area of a single element. The stabilization matrix we choose is that in [8,Sect. 3.3.2], which we briefly describe here. Consider a macroelement comprising a 2 × 2 patch of elements. Then the kth macroelement stabilization matrix is Additionally, Q k = | k |I , and the connectivity matrix that maps pressure degrees of freedom on a macroelement to global degrees of freedom is the identity, i.e. L = I . For these elements B T has full rank except in the case of Dirichlet boundary conditions, in which case null(B T ) = span{v 1 , v 2 } where v 1 and v 2 are as in (3.7). Similarly to Q 1 -Q 1 elements, we can bound Case II eigenvalues for any preconditioner that satisfies (2.8). Proof The proof is analogous to that of Corollary 3.8.

Possible extensions
It is clear that the methodology outlined above could be applied to non-square domains to ascertain the presence of Case II eigenvalues, although it may be more difficult to determine the appropriate nullspace vectors [20], and the connectivity matrix may be more complicated to describe.
Since v 2 ∈ null(B T ), −0.25/α is a Case II eigenvalue. Moreover, after boundary conditions are applied we find that (−0.75/α, v 4 ) is an additional Case II eigenvalue. For Q 1 -P 0 elements we find that (0, v 1 ), (−4/α, v 2 ), (−2/α, v 3 ) and (−2/α, v 4 ) are eigenpairs of (−C, αQ). However, because this problem has a natural outflow condition there are no Case II eigenpairs. We note that exactly the same results hold for the obstacle problem described in the next section.
Although in this section we used the ideal preconditioner P α , additional numerical experiments (described in Sect. 5), that are conducted with A replaced by a single Vcycle of algebraic multigrid (AMG), show that only λ m+1 , which takes values between 0.84 and 0.94, changes significantly when this approximation is made. Additionally, we note that we could replace diag(Q), the diagonal of the mass matrix, by lump(Q), the lumped mass matrix whose entries are the row sums of Q, or by a fixed number of iterations of Chebyshev semi-iteration. Both approaches cause λ 1 , λ p , λ m+1 and λ m+n to better approximate the values obtained when H = Q. In particular, the analysis in this section could be straightforwardly adapted for lumped mass matrices; this is particularly easy for the Q 1 pressure mass matrices considered here for which lump(Q) = 2.25 diag(Q).

Summary and interpretation
Theorem 2.1 and the subsequent analysis tells us that increasing the parameter α in P α leads to more clustered eigenvalues of P −1 α A for a range of Stokes problems, and should result in more rapid convergence of the MINRES algorithm. The likely drawback of this was that the negative eigenvalues of the preconditioned system could approach zero at a rapid rate as α is increased-in this section we have shown that this does not occur.
A key question is therefore how the value of α should be selected for practical computations. Although the theory suggests there is no "optimal" choice, we believe that a reasonable selection, and one that fits with the desire for the eigenvalues of P −1 α A to be bounded away from zero, is one that ensures that P α is as well conditioned as possible. It is well known [8, Chapter 1] that the eigenvalues of the pressure mass matrix are contained within [c m hd , C m hd ], and those of the stiffness matrix within [c a hd , C a hd −2 ], for positive h-independent constants c m , C m , c a , C a , and withd the dimension of the problem. Therefore, when seeking the best possible conditioning of P α , by "balancing" the blocks A and α H , it important to ensure that the parameter α does not exceed an O(h −2 ) value. In practice, we find that a much more moderate scaling, such as O (10), is sufficient to ensure more rapid convergence of the MINRES algorithm for a range of Stokes problems.

MINRES convergence bounds for Stokes problems
In the previous two sections we characterized the effects of α on the eigenvalues of P −1 α A. It is now of interest to ascertain the effect of varying α on the number of iterations required for preconditioned MINRES to converge to a fixed tolerance.
The following MINRES convergence bounds are well known (see, e.g., [8,Sect. 4.2.4]): where k is the set of polynomials of at most degree k and σ (P −1 Although this bound can be pessimistic, it will still provide some insight into the effect of α on preconditioned MINRES convergence.

Lemma 4.1
Let the Stokes equations be discretized by Q 2 -Q 1 elements in R 2 , and assume that (2.6) holds. Let P α be as in (1.6). Then, the eigenvalues of Alternatively, if H = diag(Q) then Here, γ is the inf-sup constant in (2. The bounds for H = diag(Q) are obtained similarly. The only additional step is in bounding ν. Since, for all y ∈ R m , y = 0, it holds that [28] Using this inequality gives the required bounds.

Lemma 4.2
Let the Stokes equations be discretized by Q 1 -Q 1 or Q 1 -P 0 elements in R 2 , and assume that (2.7) holds. Let P α be as in (1.6). Then the eigenvalues of Alternatively, for Q 1 -Q 1 elements if H = diag(Q) then, assuming that λ p = −0.25/α, Here, γ is as in ( To show that all eigenvalues are bounded below by −a, first note that since Q 1 -P 0 and Q 1 -Q 1 elements satisfy the ideal stabilization property, the largest eigenvalue of Q −1 C is less than or equal to 1. Thus, no Case II eigenvalue is smaller than −1/α. Since Case II eigenvalues are no smaller than −a. It is straightforward to show that Case III eigenvalues are also bounded below by −a. If H = diag(Q) the proof is similar to the above if we again use (4.3) to replace Q by diag(Q) in ν and μ.
To assess the effect of α on (4.2), and hence on the convergence rate of preconditioned MINRES, we compute a, b, c, d using Lemma 4.1 or 4.2 (see Tables 6 and  7). Comparison with Tables 1 and 2 shows that the eigenvalue bounds for Q 2 -Q 1 elements are very tight. For the stabilized elements a and d overestimate the magnitude of the extreme eigenvalues, but in almost all cases only by a small amount. The exception is a for Q 1 -Q 1 elements, which is close to twice |λ 1 |.
We then increase a or d so that both intervals [−a, −b] and [c, d] are of equal length, and apply (4.2). The results, in Fig. 2, clearly show that increasing α reduces η, but that as we increase α beyond about 10, η decreases much more slowly. In other words, as α is increased beyond this point, we would not anticipate a further significant reduction in iteration numbers for our preconditioned solver. This therefore motivates a value of α equal to roughly 10, as this choice essentially achieves the optimal predicted convergence rate, while at the same time ensuring that the negative eigenvalues of P −1 α A are far from zero. This pattern of behavior will be realized in numerical experiments discussed in the next section.
We end this section by discussing the effect of α on the norm used to measure convergence of preconditioned MINRES. A common stopping criterion is a specified reduction in the preconditioned residual norm, i.e. given a symmetric positive definite preconditioner P we terminate when is the kth residual and [(x In this sense increasing α relaxes the stopping criterion for the constraint equation B x = g. In our experience this is not a problem, as we show in the next section.

Numerical verification
Having motivated the application of scaled saddle-point preconditioners to Stokes problems, we would like to illustrate numerically the effect of the scaling. In particular, it is important to observe the effectiveness of this strategy when state-of-the-art preconditioners are applied both exactly and inexactly (as an inexact application will generally result in a more efficient algorithm), and determine the potency of our approach for different finite element discretizations.
To ascertain this we run the preconditioned MINRES algorithm on particular test problems, to a preconditioned residual norm tolerance of 10 −6 , within the IFISS software system [9,10,22] in matlab. In particular, we examine the regularized cavity Table 6 Eigenvalue bounds from      9 4 ] × [− 1 4 , 1 4 ] removed (see Fig. 4). No-flow conditions are applied at the top and bottom walls, and at the obstacle boundary. We impose a Poiseuille flow condition, that is v x = 1 − y 2 , v y = 0, on the inflow boundary; we also specify a natural boundary condition on the outflow boundary. In Fig. 3 we present a streamline plot for the velocity solution of the cavity problem, and a plot of the pressure solution; for these plots we set h = 2 −8 (corresponding to the finest mesh tested). In Fig. 4 we provide the same plots for the obstacle flow problem, with h = 2 −7 .  Table 8 we present iteration numbers for the MINRES solution of the regularized cavity problem using stabilized Q 1 -Q 1 elements on a uniform mesh. Within the preconditioner, we use one AMG V-cycle with point damped Gauss-Seidel smoothing for the matrix A, and 10 steps of Chebyshev semi-iteration [12,13,30] for H . We present results for different values of the (uniform) mesh parameter h, as well as values of α within P α . We observe that when α is increased, the iteration numbers clearly decrease, and there is hence a considerable benefit to applying the scaled preconditioner. This is observed for all values of mesh parameter tested. We present these results pictorially in Fig. 5, illustrating the effect of α for all values of h tested.

Effectiveness for different preconditioning options
We now wish to observe whether our approach is effective for a range of (exact and inexact) preconditioners, as well as for different finite element basis functions. In Table 10 we therefore present iteration numbers for the solution of the obstacle problem in such scenarios, with the different preconditioning strategies presented in Table 9. The matrix A is either taken to be A or an AMG V-cycle applied to it; the preconditioner H for the Schur complement is either the diagonal of Q or 10 steps of Chebyshev semi-iteration applied to Q. We highlight that we also ran the same tests with H = Q, and obtained very similar results as when using Chebyshev semiiteration. When Q 1 -P 0 finite elements are used, Q is diagonal, so we only run the tests for preconditioner options 1 and 3. In all cases the mesh parameter is fixed as h = 2 −7 , and different values of α are again taken within P α . We see that applying Chebyshev semi-iteration within the Schur complement approximation results in faster convergence than a diagonal approximation; using AMG to approximate the (1, 1)block yields roughly similar convergence as an exact inverse for Q 1 -Q 1 elements, but higher iteration counts for Q 1 -P 0 and Q 2 -Q 1 elements. Significantly, we once again observe the advantage of increasing α within the preconditioner-this behavior is replicated for all preconditioning options tested when stabilized finite elements are used. We highlight that each MINRES iteration requires the same computational operations for a given matrix system, and therefore a reduction in the iteration count results in a corresponding decrease in computing time. In the best case we observe a reduction of 30% in MINRES steps and hence CPU time, when increasing the value of α for a stabilized problem. It is important to highlight that, although the classical stopping criteria for MINRES involves convergence of the relative preconditioned residual norm to a desired tolerance, we wish to achieve accurate solutions in measures that are not themselves influenced by the preconditioner. We have therefore calculated r k 2 / r 0 2 for the solutions obtained using our solver for all values of α tested. (Recall that our stopping criterion is r k P −1 α / r 0 P −1 α < 10 −6 .) In Table 11 we state, for a range of basis functions and values of h, the 'worst case' relative residual norm achieved for the obstacle problem, and the value of α for which it was achieved. This verified that the solutions we obtained were accurate in a real sense, and not in a measure that was itself affected by the value of α. In fact, for smaller values of h (i.e. problems of higher dimension), the accuracy of our solutions seemed to improve. We observed that the rate of convergence achieved was dictated by the factor R α , as suggested by the analysis of Sect. 2.

General problems
We now investigate whether our results hold for more general problems. In particular, we examine the solution of the three-dimensional cavity problem on = [0, 1] 3 , with f = 0 and boundary conditions where v = [v x , v y , v z ] T . As for the two-dimensional cavity problem, the flow is enclosed and so the preconditioned system is singular. We also computed a two-   The G0 mesh has 22 348 degrees of freedom and the G1 mesh has 99 710 degrees of freedom dimensional Stokes flow around a circular shaped obstacle using a highly unstructured mesh, G0, (see Fig. 6) and a uniformly refinement of it, mesh G1. These results were computed using the T-IFISS software package. 5  Tables 12 and 13 show iteration numbers for different α, with the preconditioners as in Table 9. For both problems we see qualitatively similar behavior to that in Table 10 for the 2D obstacle problem, so that increasing α reduces the number of iterations needed. This is mirrored by eigenvalue computations (not shown) which, in both cases, display qualitatively similar behavior to the 2D model problems.

Concluding remarks
This work shows that including a simple scaling to well-established block diagonal preconditioners for Stokes problems can result in significantly faster convergence when applying the preconditioned MINRES method. We demonstrated theoretically why this occurs by analyzing the eigenvalues of the preconditioned matrix P −1 α A. In particular, the positive eigenvalues cluster near 1 as the scaling parameter is increased, with the negative eigenvalues also clustering and only approaching 0 slowly. We also show that the performance gains can be significant (30% reduction in CPU times) if a stabilized mixed approximation method is in use.