Constraint Schur complement preconditioners for nonsymmetric saddle point problems

Applied Numerical Mathematics 59 (2009) 151–169 www.elsevier.com/locate/apnum

Constraint Schur complement preconditioners for nonsymmetric saddle point problems ✩ Zhi-Hao Cao School of Mathematical Sciences and Laboratory of Mathematics for Nonlinear Sciences, Fudan University, Shanghai 200433, People’s Republic of China Available online 16 January 2008

Abstract We consider constraint preconditioners for block two-by-two generalized saddle point problems, this is the general nonsymmetric, nonsingular case where the (1,2) block need not equal the transposed (2,1) block and the (2,2) block may not be zero. The constraint preconditioners are derived from splittings of the (1,1) block of the generalized saddle point matrix. We show that fast convergence of the preconditioned iterative methods depends mainly on the quality of the splittings and on the effectively solving for the Schur complement systems which arise from the implementation of the constraint preconditioners. Results concerning the eigensolution distribution of the preconditioned matrix and its minimal polynomial are given. To demonstrate the effectiveness of the constraint Schur complement preconditioners we show convergence results and spectra for two model Navier–Stokes problems. © 2008 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65F10; 65N20 Keywords: Preconditioner; Saddle point problems; Schur complement; Inner-outer iterations; Navier–Stokes equations

1. Introduction We attempt to solve block 2 × 2 linear systems of the form        x f x A BT = , A ≡ C −D y g y

(1.1)

where A ∈ Rn,n is nonsingular, B, C ∈ Rm,n , m  n and D ∈ Rm,m . Recently, a large amount of work has been devoted to problems of solving large linear systems in saddle point forms which are all special cases of (1.1): in most applications we have C = B and D = 0, see, for example, [2,8– 11,13,16–19,21,24,25]; for some problems we have C = B and D = 0 but D2 is small in many such problems, the nonzero (2,2) block arises from a stabilization term, see, for example, [3,6,10,12,15,20,21,33–35,39]; finally, we note that if one uses Chebyshev collocation approach to discretize, for example, the Stokes equations, then one has a saddle point system with C = B (cf. [5,22,28,32,37]). On the solution methods for saddle point systems there is a very good survey [4]. When the coefficient matrix A above is assumed to be symmetric and D = 0, i.e. ✩

This work is supported by NSFC Projects 10171021 and 10471027. E-mail addresses: [email protected], [email protected].

0168-9274/$30.00 © 2008 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2008.01.002

152

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

 A=

A B

BT 0

 (1.2)

with A being symmetric, Keller et al. [24] investigated the use of (so called constraint) preconditioners of the form   G BT (1.3) G= B 0 with a CG-like method to solve the saddle point systems, where G is also assumed to be symmetric. They determined the eigensolution distribution of the preconditioned matrix G −1 A and an upper bound on the degree of its minimal polynomial. Cao [13] extended the constraint preconditioner to nonsymmetric saddle point matrices, i.e. A and G respectively in (1.2) and (1.3) may be nonsymmetric, while Dollar [15] extended the constraint preconditioner to symmetric saddle point matrices with (2,2) block being negative semidefinite. For the general saddle point matrix A in (1.1) Siefert and de Sturler [32] studied block diagonal preconditioners (see also [23])   G G= (1.4) D + CG−1 B T and a class of preconditioners (in general better than (1.4)) which are derived from (1.4) and referred to as related system (preconditioners) and which are of the same type as the constraint preconditioners introduced in this paper. The paper [32] is closely related to the paper [14]. In [14] a detailed analysis is provided for these two classes of preconditioners for the cases where B = C, D = 0 and A can be nonsymmetric and singular. [32] extends these preconditioners to the case where D = 0 and to allow for approximations to the Schur complement matrix that arises in the preconditioner. In this paper we study constraint preconditioners for the general saddle point systems (1.1). The constraint preconditioners are derived from splittings of the (1,1) block of A. We will show that the fast convergence of the preconditioned iterative methods with the constraint Schur complement preconditioners is mainly determined by the quality of the splittings and on the efficiency of the solution for the Schur complement systems which arise from the implementation of the constraint preconditioners. Results concerning the spectrum and the form of the eigenvectors of the preconditioned matrix G −1 A are determined. An upper bound on the degree of its minimal polynomial is given which provides an upper bound on the maximum requested number of iterations of a Krylov subspace method such as GMRES. Finally, we explore two model problems: The first, which arises from stabilized mixed finite element discretization of the Navier–Stokes equations, has D = 0 and A = AT . The second, which arises from a spectral (Chebyshev collocation) approach for an incompressible Stokes problem, has B = C. We use our theoretical results, i.e. eigenvalue bounds (cf. Theorem 2.4) and numerical experiments to illustrate that reasonable choices for splittings and effectively solving for the Schur complement systems yield good convergence. Although our eigenvalue bounds may not be sharp, they nevertheless indicate good eigenvalue clustering for reasonable choices for splittings. Note that throughout this paper  ·  indicates the 2-norm. 2. Constraint preconditioners We introduce constraint preconditioners for matrix A in (1.1) as a straightforward generalization of preconditioners in [20,24]   G BT , (2.1) G= C −D where G ∈ Rn,n is nonsingular. On the nonsingularity of A and G we have Lemma 2.1. The matrices A and G respectively in (1.1) and (2.1) are nonsingular if and only if the Schur complements −(D + CA−1 B T ) and −(D + CG−1 B T ), respectively, are nonsingular. Proof. We have

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

 A= and

 G=

I CA−1

I CG−1



A

I 

G

I

BT −(D + CA−1 B T )

153

 (2.2)

 BT , −(D + CG−1 B T )

(2.3)

hence, the conclusion of Lemma 2.1 follows immediately.

2

Remark 2.1. If D = 0, then B and C must be of full rank for A and G to be nonsingular. In contrast, if C = B, A and G are positive definite (i.e. positive real) and D is symmetric positive semidefinite, then B being of full rank implies that A and G are nonsingular. For the general discussion on the nonsingularity conditions of the saddle point matrices we refer to [4, §3.2]. In this paper we suppose that the generalized saddle point matrix A in (1.1) and the constraint preconditioner G in (2.1) satisfy all the following assumptions: • • • •

A and G are nonsingular; A and G are nonsingular or, equivalently, D + CA−1 B T and D + CG−1 B T are nonsingular; B and C are of full rank; rank(D) = p ( 0). Using (2.3) the constraint preconditioned matrix G −1 A can be expressed as   In − (In − N M)T 0 −1 G A= , −MT Im

where M = (D + CG−1 B T )−1 C, N = G−1 B T and T = I − G−1 A. Thus, G −1 A has at least m eigenvalues at 1 (cf. Theorem 2.1 below). 2.1. Eigenvalue distribution of G −1 A Now we consider the spectrum of the preconditioned matrix G −1 A. Following the approach in [2] (see also [13]) let B = U ΣV T

and C = W Γ Z T

(2.4)

be the singular value decompositions of the matrices B and C, respectively, where Σ = [Σ0 , 0] and Γ = [Γ0 , 0]

(2.5)

with Σ0 and Γ0 being m × m diagonal matrices with positive diagonal entries, U and W are orthogonal matrices of order m, and V and Z are orthogonal matrices of order n. Then we have      T  T G BT V V GZ Z ΣT , G≡ = C −D W Γ −W T DU UT     T  T  V V AZ Z ΣT A BT . = A≡ C −D W Γ −W T DU UT Let  = V T GZ, G then T ≡G

−1

 = V T AZ, A 

A=



Z U

 ΣT G  Γ −D

 = W T DU, D −1 

 ΣT A  Γ −D

(2.6) 

T

Z U

.

(2.7)

154

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

+ΓG −1 Σ T , T = I − G −1 A  = Z T (I − G−1 A)Z ≡ Z T T Z, then Let  S=D      −1  Σ T −1  −1 Σ T  G I G G S −1 = −1 I  −Γ G Γ −D − S −1 and



   −1 + G −1 Σ T   ΣT  A S −1 Γ T 0 A G  = Γ −D Im − S −1 Γ T   −1 Σ T  In − (In − G S −1 Γ )T 0 . = Im − S −1 Γ T

T ≡

 ΣT G  Γ −D

−1 

(2.8)

Theorem 2.1. The preconditioned matrix T ≡ G −1 A has at least 2m − p eigenvalues at 1. Proof. Since T ≡G

−1

 A=



Z U

 Z  T

T U

,

T and T have same eigenvalues. We know from (2.8) that T has an eigenvalue at 1 with multiplicity m and the other −1 Σ T  n eigenvalues are those of the matrix In − (In − G S −1 Γ )T.  Now we consider the eigenvalue 0 of the matrix In − G−1 Σ T  S −1 Γ . Let x be an eigenvector, then −1 Σ T  S −1 Γ )x = 0 (In − G which gives +ΓG −1 Σ T )−1 Γ x.  = Σ T (D Gx  where xˆ1 ∈ Rm , xˆ2 ∈ Rn−m , then we have Let xˆ ≡ [xˆ1T , xˆ2T ]T = Gx, +ΓG −1 Σ T )−1 Γ G −1 x, ˆ xˆ1 = Σ0 (D Let

x˜1 = Σ0−1 xˆ1 ,

xˆ2 = 0.

we have

+ΓG −1 Σ T )−1 Γ G −1 Σ T x˜1 x˜1 = (D x˜1 = 0. Since D  = W T DU, rank(D)  = rank(D) = p. Therefore, the dimension of nullspace of the which implies D −1 T −1   −1 Σ T   matrix (In − G Σ S Γ )T is at least m − p. That is to say that the matrix In − (In − G S −1 Γ )T has an eigenvalue 1 with multiplicity at least m − p. 2 −1 Σ T   ≡ In − (In − G S −1 Γ )T has at least m − p linearly From the proof of Theorem 2.1 we know that the matrix X o n,r  − In ), independent eigenvectors corresponding to the eigenvalue 1. Let T 2 ∈ R be an orthonormal basis for N (X o n,n−r   be an orthonormal basis for the complement space of N (X − In ), for the nullspace of X − In and let T 1 ∈ R o  − In )T , as T 1 . Then we have r  m − p  − In )T ), the range of (X example, we can take an orthonormal basis of R((X and    o o −1  o o  11 o −1 o X   , (2.9) T XT ≡ T 1 , T 2 X T 1 , T 2 ≡  X21 Ir 11 ∈ Rn−r,n−r , X 21 ∈ Rr,n−r . Denote Y  = − where X S −1 Γ T, we have from (2.8) that ⎞ ⎛ 11  o    o  X −1  X T T  ⎠ ⎝ I X , 21 r  Im = Y Im Im 1 Y 2 Im Y 1 ∈ Rm,n−r , Y 2 ∈ Rm,r and [Y 1 , Y 2 ] = Y T . where Y Let o

(2.10)

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

⎛

11 X   ⎝ X T = 21 1 Y

155

⎞ Ir  Y2

⎠,

(2.11)

Im

we have from (2.7), (2.8) and (2.10) that     o Z Z T T T ≡ G −1 A = U Im



−1

o

T

U

Im

(2.12)

.

We now consider the degree of the minimal polynomial of the preconditioned matrix T which determines the convergence behavior of a Krylov subspace method such as GMRES [13,24,27,29–31,38]. Lemma 2.2. Let pk (λ) be a monic polynomial of degree k and T (see (2.11)) be the matrix similar to the preconditioned matrix T ≡ G −1 A (see (2.12)), then ⎞ ⎛  11 ) (X11 − I )pk (X ⎠ 11 ) 21 pk (X (T − I )pk (T) = ⎝ X 0 k 0 k Z W k and Z k . for some matrices W Proof. If k = 1, then pk (λ) = λ − c1 . We have ⎞⎛ ⎛ 11 − c1 I X11 − I X 21 ⎠⎝ X 21 (T − I )p1 (T) = ⎝ X 0 1 2 0  Y Y Y ⎛  1 ⎞  (X11 − I )p1 (X11 ) ⎠, 11 ) 21 p1 (X =⎝ X 0   Z1 0 W1

⎞ (1 − c1 )Ir 2 Y

⎠ (1 − c1 )Im

1 = Y 1 (X 11 − c1 I ) + Y 2 X 1 = (1 − c1 )Y 2 . Thus Lemma 2.2 holds true for k = 1. 21 , Z where W Assume that Lemma 2.2 holds true for any monic polynomial with degree less than k. Any monic polynomial pk (λ) with degree k can be expressed as pk (λ) = pk−1 (λ)(λ − ck ), where pk−1 (λ) is a monic polynomial with degree k − 1. We have (T − I )pk (T) = (T − I )pk−1 (T)(T − ck I ) ⎞⎛  ⎛  11 ) (X11 − I )pk−1 (X X11 − ck I ⎠⎝ X 11 ) 21 pk−1 (X 21 =⎝ X 0 k−1 0 1 k−1 Z Y W ⎛  ⎞  (X11 − I )pk (X11 ) ⎠, 11 ) 21 pk (X =⎝ X 0 k 0 k Z W

⎞ (1 − ck )Ir 2 Y

⎠ (1 − ck )Im

k = W k−1 (X 11 − ck I ) + Z k−1 X K = (1 − ck )Z k−1 . The rest of the proof follows by induction. 21 , Z where W

2

Theorem 2.2. Let T (see (2.11)) be the block lower triangular matrix similar to the preconditioned matrix T ≡ G −1 A 11 is k (0  k  n − m + p), then the degree (see (2.12)). If the degree of the minimal polynomial of its (1,1) block X −1 of the minimal polynomial of T ≡ G A is at most k + 2 ( n − m + p + 2). 11 , then pk (X 11 ) = 0 and we have by Lemma 2.2 that Proof. Let pk (λ) be the minimal polynomial of X

0 0 0   (T − I )pk (T ) = 0 0 0 , k Z k 0 W

156

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

thus (T − I )2 pk (T) = 0. Since T and T are similar (cf. (2.12)), they have same minimal polynomial. The theorem follows.

2

Remark 2.2. If in (2.11) r = 0 which implies rank(D) ≡ p = m, i.e. D is of full rank, T in (2.11) is simplified to   11 X  . T =  Y1 Im 11 , then In this case if k is the degree of the minimal polynomial pk (λ) of X   11 ) 0 11 − I )pk (X (X = 0. (T − I )pk (T) = 11 ) 1 pk (X 0 Y Therefore, the degree of the minimal polynomial of T ≡ G −1 A is at most k + 1. This result is consistent with the main theoretical results, i.e. Propositions 2.1, 2.2 and 3.1 in [1]. From (2.4) and (2.5) we have BV ≡ B[V1 , V2 ] = [U Σ0 , 0], CZ ≡ C[Z1 , Z2 ] = [W Γ0 , 0], i.e. V2 and Z2 are respectively orthonormal bases of N (B) and N (C). Let  T   P1 p T = Q1 p P1T D = QP ≡ [Q1 , Q2 ] 0 P2T

(2.13)

(2.14)

be the singular value decomposition of the matrix D, where Q ≡ [Q1 , Q2 ] and P ≡ [P1 , P2 ] are orthogonal matrices of order m, Q1 , P1 ∈ Rm,p and Q2 , P2 ∈ Rm,m−p , while p is a p × p diagonal matrix with positive diagonal entries. We note that D T Q1 = P1 p ,

DP2 = 0,

i.e. N (D) = R(P2 ), R(D T ) = R(P1 ). From (2.6) we have ⎛ ⎞  T      V1T AZ1 V1T AZ2 Σ0 T A B V ⎝ V T AZ1 V T AZ2 ⎠ Z A≡ , = 0 2 2 C −D W UT T   0 −Q1 p P1 Γ0 ⎞ ⎛ T T GZ  T      GZ V Σ0 V 1 2 1 1 T Z G B V ⎠ ⎝ V T GZ1 V T GZ2 G≡ , = 0 2 2 C −D W UT T   0 −Q1 p P1 Γ0

(2.15)

(2.16)

where 1 = W T Q1 , Q

1 = U T P1 . P

(2.17)

We now consider the eigenvector distribution of the preconditioned matrix T ≡ G −1 A. We have the following theorem which reduces to Theorem 4 in [13] when the (2,2) block in A and in G is zero (i.e. D = 0) and C = B. Theorem 2.3. The preconditioned matrix G −1 A has m + i + j + k linearly independent eigenvectors. There are (i) m + i (0  i  n) eigenvectors corresponding to the eigenvalue 1 (cf. (2.23)), where i = dimN (A − G); (ii) j (0  j  n − m) eigenvectors corresponding to eigenvalues λl , l = 1, . . . , j, not being 1 (cf. (2.26)), where j is the number of linearly independent eigenvectors of (V2T GZ2 )−1 (V2T AZ2 ) corresponding to the eigenvalues λl , l = 1, . . . , j , where V2 and Z2 are orthonormal bases of N (B) and N (C), respectively;

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

157

(iii) k (0  k  p) eigenvectors corresponding to eigenvalues not being 1 (cf. (2.31)), where k is the number of T −1 T −1 T linearly independent eigenvectors of the matrix (G + B T P1 −1 p Q1 C) (A + B P1 p Q1 C) corresponding to eigenvalues not being 1, where P1 , p and Q1 are defined in (2.14). Proof. Let us consider the eigenvalue problem −1         A BT x x x G BT =λ . G −1 A ≡ C −D C −D y y y By (2.16) this eigenvalue problem is equivalent to the generalized eigenvalue problem: ⎛ T ⎞ ⎞ ⎛ T V1 AZ1 V1T AZ2 Σ0 V1 GZ1 V1T GZ2 Σ0 xˆ1 xˆ1 ⎝ V T AZ1 V T AZ2 ⎠ xˆ2 = λ ⎝ V T GZ1 V T GZ2 ⎠ xˆ2 , 0 0 2 2 2 2 p p yˆ yˆ 0 − 0 − Γ0 Γ0

(2.18)

1 p P T . We have from (2.18) p = Q where xˆ1 = Z1T x, xˆ2 = Z2T x, yˆ = U T y and  1 ˆ V1T AZ1 xˆ1 + V1T AZ2 xˆ2 + Σ0 yˆ = λ(V1T GZ1 xˆ1 + V1T GZ2 xˆ2 + Σ0 y),

(2.19)

V2T AZ1 xˆ1

(2.20)

+ V2T AZ2 xˆ2 = λ(V2T GZ1 xˆ1 + V2T GZ2 xˆ2 ), 1 p P 1 p P 1T yˆ = λ(Γ0 xˆ1 − Q 1T y). Γ0 xˆ1 − Q ˆ

(2.21) T

Obviously, there are m linearly independent eigenvectors [0T , 0T , yˆ (l) ]T , l = 1, . . . , m, corresponding to the eigenvalue 1. For λ = 1 (2.21) holds trivially while (2.19) and (2.20) are simplified to ˆ V T AZ xˆ = V T GZ x,

(2.22)

where xˆ = [xˆ1T , xˆ2T ]T . (2.22) is equivalent to xˆ ∈ N ((A − G)Z). Let i = dim N ((A − G)Z) (0  i  n). If i > 0, let (l)T

(l)T

xˆ (l) ≡ [xˆ1 , xˆ2 ]T ∈ N ((A − G)Z), l = 1, . . . , i, be the linearly independent nullvectors of (A − G)Z. Then each of the vectors  (l)T (l)T T ∈ Rn+m , l = 1, . . . , i, xˆ , zˆ where zˆ (l) ∈ Rm is arbitrary, is an eigenvector of T corresponding to the eigenvalue 1. It is easy to see that there are m + i linearly independent eigenvectors of T corresponding to the eigenvalue 1 which can be taken as  T T (l)T T  (l)T (l)T (l)T T 0 , 0 , yˆ , l = 1, . . . , m; xˆ1 , xˆ2 , zˆ , l = 1, . . . , i, (2.23) where yˆ (l) ∈ Rm , l = 1, . . . , m, are arbitrary m linearly independent vectors, zˆ (l) ∈ Rm , l = 1, . . . , i, are arbitrary i vectors. Note that dim N ((A − G)Z) = dim N (A − G). If λ = 1, then from (2.21) we have 1 p P 1T y. ˆ Γ0 xˆ1 = Q

(2.24)

T yˆ ≡ P T U yˆ = P T y = 0, i.e. y = P2 P T y ∈ N (D) (cf. (2.15)), then (2.24) implies xˆ1 = 0, (2.19) and Case 1. If P 1 1 1 2 (2.20) are simplified to ˆ V1T AZ2 xˆ2 + Σ0 yˆ = λ(V1T GZ2 xˆ2 + Σ0 y), V2T AZ2 xˆ2 = λV2T GZ2 xˆ2 .

(2.25)

It is easy to see from (2.25) that if xˆ2 = 0, then yˆ = 0. Thus xˆ2 = 0. Suppose is nonsingular, then from (2.25) we know that if [0T , xˆ2T , yˆ T ]T is an eigenvector of T corresponding to an eigenvalue λ not being 1, then the vector xˆ2 is 1  T −1 T an eigenvector of (V2T GZ2 )−1 (V2T AZ2 ) corresponding to the eigenvalue λ and yˆ = λ−1 P2 P2 Σ0 V1 (A − λG)Z2 xˆ2 T T T 2 = U P2 ), since y = P2 P y which is equivalent to yˆ = P 2 P  y. (where P ˆ Let j (0  j  n − m) be the number 2 2 T −1 of the linearly independent eigenvectors of (V2 GZ2 ) (V2 AZ2 ) corresponding to eigenvalues not being 1, then j linearly independent eigenvectors of T can be taken as V2T GZ2

158

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

 T (l)T (l)T T 0 , wˆ 2 , qˆ ,

l = 1, . . . , j,

(2.26)

where wˆ 2 , l = 1, . . . , j , are linearly independent eigenvectors of the matrix (V2T GZ2 )−1 (V2T AZ2 ) corresponding to eigenvalues λl = 1, l = 1, . . . , j , and (l)

1  T −1 T (l) P2 P2 Σ0 V1 (A − λl G)Z2 wˆ 2 , l = 1, . . . , j. λl − 1 T yˆ ≡ P T y = 0, i.e. y = P1 y1 + P2 y2 with y1 = 0, in this case xˆ1 = 0 (cf. (2.24)). Let y = P1 y1 , i.e. Case 2. If P 1 1 1 P T yˆ = y. ˆ Then (2.24) implies y ∈ R(D T ) (cf.(2.15)), we have (cf.(2.17)) P qˆ (l) =

1

yˆ

T 1 −1 =P p Q1 Γ0 xˆ 1 .

(2.27)

Substituting (2.27) into (2.19) gives T T 1 −1 T  −1 T V1T AZ1 xˆ1 + V1T AZ2 xˆ2 + Σ0 P p Q1 Γ0 xˆ 1 = λ(V1 GZ1 xˆ 1 + V1 GZ2 xˆ 2 + Σ0 P1 p Q1 Γ0 xˆ 1 ).

(2.28)

Noting that (cf. (2.17) and (2.13)) T −1 T T T −1 T 1 −1 T Σ0 P p Q1 Γ0 xˆ 1 = Σ0 U P1 p Q1 W Γ0 xˆ 1 = V1 B P1 p Q1 CZ1 xˆ 1 ,

(2.28) can be rewritten as

  T T T −1 T V1T (A + B T P1 −1 p Q1 C)Z1 xˆ 1 + V1 AZ2 xˆ 2 = λ V1 (G + B P1 p Q1 C)Z1 xˆ 1 + V1 GZ2 xˆ 2 .

Combining which with (2.20) gives      T T T xˆ1 V1 B P1 −1 p Q1 CZ1 0 V T AZ + xˆ2 0 0      T T T −1 xˆ1 V1 B P1 p Q1 CZ1 0 . = λ V T GZ + xˆ2 0 0

(2.29)

Since V2T B T = 0 and CZ2 = 0 (cf.(2.13)), we have  T T  T V1 B P1 −1 p Q1 CZ1 0 = V T B T P −1 QT CZ. 1 p 1 0 0 Thus (2.29) can be expressed as follows. T T T −1 T V T (A + B T P1 −1 p Q1 C)Z xˆ = λV (G + B P1 p Q1 C)Z xˆ

which is equivalent to T T −1 T (A + B T P1 −1 p Q1 C)x = λ(G + B P1 p Q1 C)x.

(2.30)

T Suppose that G + B T P1 −1 p Q1 C is nonsingular, the generalized eigenvalue problem (2.30) defines at most p linearly independent eigenvectors corresponding to eigenvalues not being 1, since (2.30) comes from y = P1 P1T y (cf. (2.27) and (2.19)) or, equivalently, from y ∈ R(D T ) (cf. (2.15)). Let k (0  k  p) be the number of linearly independent eigenvectors of (2.30) corresponding to eigenvalues not being 1, then k linearly independent eigenvectors of T can be taken as  (l)T (l)T (l)T T uˆ 1 , uˆ 2 , vˆ , l = 1, . . . , k, (2.31) (l)

(l)

where uˆ 1 = Z1T u(l) , uˆ 2 = Z2T u(l) , l = 1, . . . , k, while u(l) , l = 1, . . . , k, are k linearly independent eigenvectors of (2.30) corresponding to eigenvalues not being 1 and (cf. (2.27)) 1 −1 T ˆ (l) , vˆ (l) = P p Q1 Γ0 u 1

l = 1, . . . , k,

are k linearly independent vectors. 1 ), l = 1, . . . , k, are linearly independent and qˆ (l) ∈ R(P 2 ), l = 1, . . . , j , the k eigenvectors Since vˆ (l) ∈ R(P in (2.31) and the j eigenvectors in (2.26) are linearly independent. The m + i eigenvectors in (2.23) are corresponding to the eigenvalue 1. Thus, all the m + i + j + k eigenvectors of T are linearly independent.

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

By (2.7) and (2.8) we have    Z Z T ≡ G −1 A = T U

T U

,

therefore, T and T have same eigenvalues. If ditioned matrix T

≡ G −1 A.

159

 xˆ 

2

yˆ

is an eigenvector of T, then

 Z xˆ  U yˆ

is an eigenvector of the precon-

2.2. Perturbation bounds on the eigenvalues of G −1 A We know from (2.7) and (2.8) that each eigenvalue λ of T ≡ G −1 A satisfies   −1 Σ T  −1 Σ T  S −1 Γ )T  In − G S −1 Γ T |λ − 1|  (In − G −1 Σ T  = In − G S −1 Γ T .

(2.32)

If G is a good approximation of A, T  ≡ In − G−1 A will be small. In this case G −1 A has all eigenvalues clustered −1 Σ T  around 1 if In − G S −1 Γ  is not large. −1 Σ T  In order to obtain an upper bounds of In − G S −1 Γ  we proceed similarly to [14,32]. Without loss of  =  ∈ Rn,m and M  ∈ Rm,n . Assume that =G −1 Σ T and M S −1 Γ , we have N generality we assume n − m  m. Let N −1 Σ T ≡ (D +ΓG −1 Σ T )−1 Γ G −1 Σ T has full set of eigenvectors N ≡ the matrix M S −1 Γ G N vˆj = (1 + θj )vˆj , M

j = 1, . . . , m.

(2.33)

 is of full rank, N V   = [vˆ1 , . . . , vˆm ] ∈ C m,m is the eigenvector matrix with vˆj  = 1, j = 1, . . . , m. Note that N Let V H m,m          is a is also a matrix of full rank. Let N V = U2 R be the QR decomposition of N V , i.e. U2 U2 = Im , R ∈ C nonsingular upper triangular matrix. (2.33) can be expressed in a matrix form  N) V =V (Im + Θ),  (M

(2.34)

 = diag(θ1 , . . . , θm ). where Θ   1 = 0 , we have (N M)  U 1 = 0. It is easy to see from (2.34) that Let U In−m M N )V R −1 = N V (I + Θ) −1 = U 2 R(I  + Θ) −1 . M)  U 2 = N( R R (N Thus we have 2 ] = [U 1 , −U 2 R 1 , U 2 ] Θ −1 ] = [U M][  U 1 , U R [I − N



In−m

(2.35) 

Θ −1 R −R

.

1 , U 2 ] ∈ C n,n is nonsingular, (2.36) implies Assume that [U   In−m   −1     I − N M = [U1 , U2 ] Θ −1 [U1 , U2 ] . R −R

(2.36)

(2.37)

M  we need the singular value decomposition of the matrix U H U 2 (cf. [14,32]) To analyze I − N 1 2 = Φ 1H U Ω Ψ H ≡ [φˆ 1 , . . . , φˆ n−m ] diag(ωˆ 1 , . . . , ωˆ m )[ψˆ 1 , . . . , ψˆ m ]H , U

(2.38)

1 , U 2 ] is assumed to be nonsingular. We have from (2.38) that where 1 > ωˆ 1  · · ·  ωˆ m  0, since [U 1 Φ 2 Ψ 1H (U 1 U ) = U Ω  U which implies 1 Φ 1 U 2 Ψ 2 Ψ 1H )U =U Ω  + (I − U . U Let 1 U 2 Ψ = W  ≡ [w˜ 1 , . . . , w˜ m ], 1H )U (I − U

(2.39)

160

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

H W  = 0 and from (2.39) we have then U 1 2 ψˆ j = ωˆ j U 1 φˆ j + w˜ j , U

j = 1, . . . , m,

(2.40)

which implies that ωˆ j2 + w˜ j 2 = 1, j = 1, . . . , m. Denote σˆ j = w˜ j , j = 1, . . . , m and let wˆ j = wˆ jH wˆ k = 0,

w˜ j σˆ j

, j = 1, . . . , m. Then using (2.40) and (2.38) it is easy to see that

k = j, k, j = 1, . . . , m.

Thus, we can express (2.40) in the following matrix form 2 Ψ 1 Φ Σ =U Ω +W  U

(2.41)

 ≡ [wˆ 1 , . . . , wˆ m ] ∈ C n,m , W H W H W  = Im , U  = 0.  = diag(σˆ 1 , . . . , σˆ m ), ωˆ 2 + σˆ 2 = 1, j = 1, . . . , m, and W where Σ j j 1 2 can be expressed as By (2.41) U    Ω H      H      U2 = (U1 Φ Ω + W Σ)Ψ = ( U1 Φ W )  Ψ Σ 1 , U 2 ] from which we easily obtain an expression of the matrix [U  H     Φ 1 , U  ] In−m Ω 2 ] = [U 1 Φ,  W [U  H . Σ Ψ 1 Φ,  ] is an unitary matrix, (2.42) implies  W Since [U   H T    U Σ −1  Φ In−m −Ω Φ −1 1   . [U1 , U2 ] =  −1 H Ψ Σ W 2 = I − Ω  it is easy to see that T Ω, Using (2.42), (2.43) and noting that Σ         In−m [U   1 , U 2 ]−1  = [In−m , −Ω Σ −1 ] [U  1 , U2 ]  0 1  1  −1 Ω  2 = (1 − ωˆ 12 )− 2 .  m−Ω T Ω) = In−m + Ω(I   1 , U 1 , U 2 ] 0  −1 [U 2 ]−1  as follows Following [32] we now establish a bound on [U −R Θ R       1 , U 1 , U 2 ] 0  −1 [U 2 ]−1 x [U   0 − R Θ R −1 [U  1 , U  2 ]  = max [ U , U ] 1 2   x=0 Θ −1 R −R x   −1 x2    U2 R Θ R = max . 1 x1 + U 2 x2  1 x1 +U 2 x2 =0 U U

(2.42)

(2.43)

(2.44)

(2.45)

2 R Θ −1 . We have from (2.41) that Θ −1 x2   R R R Obviously, x2 = 0, we may assume x2  = 1, so that U   1 x1 + U 2 x2  =  Σ Ω Ψ  H x2 ) + W Ψ  H x2  U U1 (x1 + Φ Ω Ψ H x2 . This gives which is minimized for any x2 by x1 = −Φ 1 x1 + U 2 x2  = W Σ Ψ  H x2  = Σ Ψ  H x2  U which is minimized for x2 = ψˆ 1 . Therefore, we have from (2.45) that       −1  2 − 1   −1 [U 1 , U  2 ] 0  [ U , U ] 1 2 −1  (1 − ωˆ 1 ) 2 R Θ R .     −R Θ R Thus, we have established a perturbation bound on the eigenvalues of G −1 A.

(2.46)

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

161

Theorem 2.4. Let T (see (2.8)) be the 2 × 2 block lower triangular matrix similar to the preconditioned matrix T ≡ −1 Σ T  =G −1 Σ T , M = S −1 Γ . G −1 A (see (2.7) and (2.8)), its (1,1) block is In − (In − G S −1 Γ )T (see (2.8)). Let N −1 Σ T has a full set of eigenvectors V  ∈ C m,m and suppose [U 1 , U 2 ] is nonsingular, N ≡ Suppose that M S −1 Γ G   1 = 0 , U 2 = N V R −1 , U H U  is an m × m nonsingular upper triangular matrix, i.e. U 2 R 2 = Im , R  is a where U 2 In−m −1 V , then each eigenvalue λ of the preconditioned matrix G A satisfies |λ − 1|  QR decomposition of the matrix N − 12 2 −1    (1 − ωˆ 1 ) (1 + R Θ R )T , where Θ is defined in (2.34) and ωˆ 1 is defined in (2.38). 3. Implementation There are various strategies that can be used to implement the proposed preconditioner in (2.1), i.e. to choose G and to compute G −1 v for a given vector v. Since G in (2.1) is assumed to be nonsingular, a block LU factorization of G can be written as    G I G−1 B T (3.1) G= I C −S where S = D + CG−1 B T . If G is taken as an incomplete LU factorization of A [26,29] A = LU + R,

G = LU,

(3.2)

where L and U are respectively lower and upper triangular matrices. Then z = G −1 v or, equivalently, the solution z = [z1T , z2T ]T of the linear system     v1 z (3.3) G 1 = z2 v2 can be obtained by the following algorithm: Algorithm 3.1. For a given vector v = [v1T , v2T ]T compute the vector z = [z1T , z2T ]T = G −1 v (1) (2) (3) (4) (5)

z1 = U −1 (L−1 v1 ); t2 = Cz1 − v2 ; Solve (D + CU −1 L−1 B T )z2 = t2 ; t1 = U −1 (L−1 (B T z2 )); z1 = z1 − t1 .

We will compare our constraint Schur complement preconditioner with the block triangular Schur complement preconditioner   G BT  (3.4) G= −(D + CG−1 B T ) (cf. [19,27]) and the block diagonal Schur complement preconditioner   G  G= D + CG−1 B T

(3.5)

presented in [32]. It should be pointed out that [32] also proposes a constraint (and approximate constraint) precondi−1 v and z = G −1 v can be obtained by the following algorithms. tioner. Using (3.2) z = G −1 v Algorithm 3.2. For a given vector v = [v1T , v2T ]T compute the vector z = [z1T , z2T ]T = G (1) Solve (D + CU −1 L−1 B T )z2 = −v2 ; (2) z1 = U −1 (L−1 (v1 − B T z2 )); −1 v Algorithm 3.3. For a given vector v = [v1T , v2T ]T compute the vector z = [z1T , z2T ]T = G

162

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

(1) z1 = U −1 (L−1 v1 ); (2) Solve (D + CU −1 L−1 B T )z2 = v2 . From these algorithms we can see that the key to obtaining fast convergence for the preconditioned GMRES lies with the choices of LU (i.e. the quality of the splitting of the (1,1) block A in (3.2)) and with the effective solution for the Schur complement system (D + CU −1 L−1 B T )y = w.

(3.6)

We can use, for example, GMRES as an inner iterative method to solve the Schur complement system (3.6). The resulting preconditioned iterative method may attain the required iteration precision by small outer iteration number, but the total inner iteration number may be very big. Thus, the whole preconditioned iterative method is timing cost. In the numerical experiment in the next section, we made use of a relaxed inner iteration stopping criteria (cf. [7,36]) to reduce the computational cost: the inner iteration residual rˆj = w − (D + CU −1 L−1 B T )y (ν) (ν)

in the j -th outer iteration satisfies (y (0) = 0) (ν)

ˆrj  w



η min(1, r (j −1) )

(3.7)

for some prescribed tolerance η, where r (j −1) is the (j − 1)-th outer iteration residual. 4. Numerical experiments As in [32, §5] we present numerical experiments for two model problems, both arising from the Navier–Stokes equations. The first model problem involves a stabilized mixed finite element discretization of the Navier–Stokes equations. The resulting saddle point system is nonsymmetric but has B = C and A being positive real. As Siefert and de Sturler pointed out in [32] that excellent work has been done by others on preconditioners for this specific problem [17,19, 35,39], which we do not intend to supplant. Our goal is to illustrate the effect of the constraint Schur complement preconditioners proposed in this paper on the convergence behavior and the eigenvalue distributions for a problem which is well-understood. The second model problem which was suggested to the author by E. de Sturler involves a spectral collocation discretization for the incompressible Stokes equations on a square [5,22,28,32,37]. In [32] this problem is also studied for a constraint preconditioner and a block diagonal preconditioner. This application has B = C and A being nonsymmetric and not positive real. As in [32] for this application, we present preconditioned GMRES convergence results as well as the locations of the eigenvalues of the preconditioned matrix. 4.1. Navier–Stokes with stabilized mixed finite element discretization The problem under consideration arises from the linearization of the steady-state Navier–Stokes equations, i.e. the Oseen problems of the following form −νu + w · u + grad p = f, div u = 0 in Ω

(4.1)

with suitable boundary conditions on ∂Ω, where Ω ⊂ R2 is a bounded domain and w is a given divergence free field. The scalar ν is the viscosity, the vector field u represents the velocity, and p denotes the pressure. The test problem is a “leaky” two-dimensional lid-driven cavity problem in a square domain (0  x  1 : 0  y  1). The boundary conditions are ux = uy = 0 on the three fixed walls (x = 0, y = 0, x = 1), and ux = 1, uy = 0 on the moving wall (y = 1). On the constructing coefficient matrix A, we take the constant wind: wx = 1, wy = 2 and the circulating wind: wx = 8x(x − 1)(1 − 2y), wy = 8(2x − 1)y(y − 1) (cf. [19]).

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

163

To discretize (4.1), we take a finite element subdivision based on ne × ne uniform grids of square elements. the mixed finite element used is the bilinear-constant velocity-pressure: Q1 − P0 pair with global stabilization or local stabilization (cf. [33,34]). If we take node bases, then the coefficient matrix A of the linear system, which is equivalent to (1.1) (with C = B), is the following ⎞ ⎛ F B1T (4.2) A=⎝ F B2T ⎠ , B1 B2 −D  + N . Then it is easy to see for the constant wind that the element matrices are following (denote where F = ν A h = 1/ne) (cf. [10,19,21]) ⎛ ⎞ ⎛ ⎞ 4 −1 −1 −2 0 2 4 3 h ⎜ −2 0 −1 4 −2 −1 ⎟ 1 4 ⎟ e = 1 ⎜ Ne = A ⎝ ⎠, ⎝ ⎠, −1 −2 4 −1 −4 −1 0 2 6 12 −2 −1 −1 4 −3 −4 −2 0 h h (B1 )e = ( 1 −1 1 −1 ) , (B2 )e = ( 1 1 −1 −1 ) (4.3) 2 2 here, the order of the vertices in the unit square element is from the left to the right and from below to above. Matrix D in (4.2) results from the stabilization via the global jump formulation D(uh , p h ) = βh

Ns  

[[uh ]][[p h ]] ds,

(4.4)

i=1∂Ω i

where β(> 0) is a stabilizing parameter, [[·]] is the jump operator, the summation is over all interior inter-element edges {∂Ωi | i = 1, . . . , Ns }, or from the stabilization via local jump formulation 4  m   D(uh , p h ) = βh [[uh ]][[p h ]] ds, (4.5) j =1 i=1∂Ω ji

where the first summation is over all macroelements of 2 × 2 elements (we assume ne is a multiple of four, hence the elements in our triangulization can be assembled into Nm disjoint macroelements) and the second summation runs over all inter-element edges strictly within each macroelement {∂Ωmi | i = 1, . . . , 4}. For a macroelement the local stabilization matrix is the following: ⎛ ⎞ 2 −1 −1 0 0 −1 ⎟ ⎜ −1 2 De = βh2 ⎝ (4.6) ⎠. −1 0 2 −1 0 −1 −1 2 When ne = 3 (h = 1/3), matrix D of the global stabilization is the following: ⎞ ⎛ 2 −1 −1 −1 ⎟ ⎜ −1 3 −1 ⎟ ⎜ −1 2 −1 ⎟ ⎜ ⎟ ⎜ 3 −1 −1 ⎟ ⎜ −1 ⎟ 2⎜ (4.7) D = βh ⎜ −1 −1 4 −1 −1 ⎟. ⎟ ⎜ −1 −1 3 −1 ⎟ ⎜ ⎟ ⎜ −1 2 −1 ⎟ ⎜ ⎠ ⎝ −1 −1 3 −1 −1 −1 2 In our numerical experiments we take β = 1 and β = 0.25 for global and local respectively (cf.  stabilizations,  [33,34]). We note that in this method of discretization, in the matrix A (2.1) A = F F , C = B = [B1 , B2 ] and D is symmetric positive semi-definite.

164

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

Table 1 Values of m and n, nonzeros in A, B, D1 and D2, order of A h 1 16 1 32 1 64

order of A

nz(A)

nz(B)

nz(D1)

nz(D2)

256

450

3698

3600

1216

768

706

1024

1922

16562

15376

4992

3072

2946

4096

7938

69938

63504

20224

12288

12034

m

n

All the runs were done in MATLAB on an IBM PC with a 2.13 GHz processor. Table 2 Iterations for circulating wind, ν = 0.1 global τ 0.1

0.01

0.005

0.001 10−4 10−5

local

1 16

1 32

1 64

1 16

1 32

1 64

13 [445] (1.22) 5 [134] (0.53) 4 [110] (0.45) 3 [91] (0.44) 2 [76] (0.33) 2 [72] (0.30)

19 [735] (6.58) 9 [256] (2.58) 7 [193] (1.98) 4 [116] (1.38) 3 [93] (1.42) 2 [79] (1.14)

34 [1711] (60.92) 16 [620] (23.6) 12 [394] (16.31) 7 [196] (10.72) 3 [101] (8.61) 3 [91] (9.92)

13 [463] (1.23) 5 [148] (0.58) 4 [128] (0.44) 3 [108] (0.45) 2 [89] (0.38) 2 [84] (0.41)

19 [709] (6.31) 9 [257] (2.59) 7 [200] (2.03) 4 [133] (1.67) 3 [108] (1.53) 2 [91] (1.16)

33 [1661] (57.81) 16 [594] (22.66) 12 [386] (15.63) 7 [200] (11.08) 3 [117] (9.91) 3 [103] (10.98)

We test the constraint Schur complement preconditioner (2.1) (with C = B), the block triangular Schur complement preconditioner (3.4) (with C = B) and the block diagonal Schur complement preconditioner (3.5) (with C = B) with GMRES. We use left preconditioning. The stopping criterion is r (k)   10−6 , r (0) 

(4.8)

the initial guess is the vector with all ones. The triangular matrices L and U in Algorithms 3.1–3.3 are obtained from an incomplete LU factorization of A with the drop tolerance [29] τ = 0.01. The inner iteration tolerance η in the relaxed stopping criteria (3.7) for solving the Schur complement system is taken as η = 10−8 . We take three values of 1 1 1 , 32 , 64 which are corresponding to three dimensions; (m = 16 × 16, n = 2 × 15 × 15), (m = 32 × 32, n = h: h = 16 2 × 31 × 31), (m = 64 × 64, n = 2 × 63 × 63). The orders of the matrix A and the number of nonzeros in the associated matrices A, B, D1 (corresponding to the global stabilization) and D2 (corresponding to the local stabilization) are given in Table 1. In the first test, following Theorem 2.4 we consider the relationship between the convergence behavior and the quality of splittings of the (1,1) block A. Therefore, we take six values of drop tolerance τ : τ = 0.1, 0.01, 0.005, 0.001, 10−4 and 10−5 to perform six incomplete LU factorizations of A, i.e. we take six splittings of the (1,1) block A. In Table 2 we show the outer iteration counts and the total inner iteration counts (each total inner iteration number is in a pair of brackets [ ] under the corresponding outer iteration number) of the constraint Schur complement preconditioned GMRES for the circulating wind with ν = 0.1. The left is for global stabilization and the right is for local stabilization. The CPU time (in second) is also showed in a pair of brackets ( ) under the total inner iteration number. From Table 2 we can see

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

165

Fig. 1. Eigenvalues of G −1 A, circulating wind, global Stabilization, ν = 0.1.

Table 3 Iterations for global stabilization, circulating wind, τ = 0.01 Const.

B-triang.

B-diag.

h

1/16

1/32

1/64

1/16

1/32

1/64

1/16

1/32

1/64

ν=1 ν = 1/10 ν = 1/20 ν = 1/30 ν = 1/50 ν = 1/100

7 5 5 5 5 5

10 9 8 8 8 7

16 16 17 16 15 14

13 13 13 13 12 12

20 21 20 20 19 19

31 34 33 33 33 34

21 22 22 22 22 22

28 31 29 28 28 28

44 52 53 51 48 45

• The iteration counts decrease with τ , i.e. the better the splitting of the (1,1) block A, the fast the convergence of the preconditioned GMRES. • The CPU time may not always decrease with τ , since the fill in of the incomplete LU factorization is increasing with τ decreasing. Thus, one should choose τ reasonably. Fig. 1 plots the eigenvalues of the preconditioned matrix G −1 A arising from the global stabilization for the circulating wind and ν = 0.1. The four subplots are corresponding to τ = 0.1, 0.01, 0.005 and 0.001, respectively. From Fig. 1 we can see • For each τ the eigenvalues cluster around the point 1. • The better the splitting of the (1,1) block A, the more tightly the eigenvalues of the preconditioned matrix cluster around the point 1. 1 1 1 1 1 In the second test, for each h we test six values of ν: ν = 1, 10 , 20 , 30 , 50 , 100 . We note that the mesh Reynolds hw number is defined as Rey = ν . Thus, for constant wind, w = 2, we have six values of the mesh Reynolds number Rey = 2h, 20h, 40h, 60h, 100h and 200h. Tables 3–6 show the (outer) iteration counts of GMRES with the constraint Schur complement preconditioner (Algorithm 3.1), with the block triangular Schur complement preconditioner (Algorithm 3.2) and with the block diagonal Schur complement preconditioner (Algorithm 3.3). Among them Tables 3 and 4 are both for the global

166

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

Table 4 Iterations for global stabilization, constant wind Const.

B-triang.

B-diag.

h

1/16

1/32

1/64

1/16

1/32

1/64

1/16

1/32

1/64

ν=1 ν = 1/10 ν = 1/20 ν = 1/30 ν = 1/50 ν = 1/100

6 5 5 4 4 4

10 8 7 6 5 4

16 12 11 10 9 6

13 11 9 8 8 9

20 18 15 14 10 9

31 27 24 22 20 11

21 24 26 24 22 22

28 28 25 26 24 21

42 41 38 35 35 24

Table 5 Iterations for local stabilization, circulating wind Const.

B-triang.

B-diag.

h

1/16

1/32

1/64

1/16

1/32

1/64

1/16

1/32

1/64

ν=1 ν = 1/10 ν = 1/20 ν = 1/30 ν = 1/50 ν = 1/100

6 5 5 5 5 5

10 9 8 8 8 7

16 16 17 16 15 14

14 14 13 13 12 12

21 21 20 20 19 19

31 34 33 33 33 32

24 22 20 20 18 18

32 30 29 28 26 27

50 54 53 51 48 42

Table 6 Iterations for local stabilization, constant wind Const.

B-triang.

B-diag.

h

1/16

1/32

1/64

1/16

1/32

1/64

1/16

1/32

1/64

ν=1 ν = 1/10 ν = 1/20 ν = 1/30 ν = 1/50 ν = 1/100

6 5 5 4 4 4

10 8 7 6 5 4

15 12 11 10 9 6

15 11 10 8 8 9

21 18 15 14 10 9

32 27 24 23 20 11

24 22 20 19 16 17

33 28 26 29 20 18

49 41 38 35 34 21

stabilization corresponding respectively to the circulating wind and to the constant wind. Analogously, Tables 5 and 6 are both for the local stabilization corresponding respectively to the circulating wind and to the constant wind. Since the main computational work in all these three methods is the same, i.e., the work of the inner iteration for solving the Schur complement system (D + B(LU )−1 B T )y = w, the (outer) iteration counts represent the computational performance of the preconditioned GMRES. From these tables we can see that • The convergence behavior for the GMRES method with the constraint Schur complement preconditioner is better than that with the block triangular Schur complement preconditioner and that with the block diagonal Schur complement preconditioner. • (Outer) iteration counts for all the preconditioners are almost independent of the viscosity ν. Indeed, in most cases they actually decrease slightly with ν. However, it should be pointed out that the total inner iteration counts often increase with ν decreasing. Thus, the whole convergence performance of the preconditioned GMRES gets worse with ν decreasing. • (Outer) iteration counts for all the preconditioners increase slowly with h−1 . 4.2. Incompressible Stokes with Chebyshev collocation We consider the homogeneous Stokes boundary value problem

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

167

Table 7 Iterations for constraint Schur complement preconditioned GMRES τ 0.1

0.01

0.005

0.001 10−4

10(319)

15(704)

20(1239)

25(1924)

30(2759)

8 [451] (0.91) 4 [247] (0.64) 4 [245] (0.64) 3 [190] (0.55) 2 [147] (0.53)

11 [892] (2.52) 4 [347] (1.56) 4 [365] (1.56) 3 [293] (1.78) 2 [222] (2.08)

10 [1064] (8.88) 5 [568] (6.39) 4 [474] (6.05) 3 [378] (6.99) 2 [288] (8.50)

12 [1582] (23.69) 6 [830] (16.36) 5 [691] (15.64) 3 [473] (15.36) 2 [357] (19.53)

13 [2031] (48.97) 8 [1228] (38.45) 6 [948] (33.75) 3 [548] (28.44) 2 [411] (37.47)

Fig. 2. Eigenvalues of G −1 A, N = 15, order = 704.

−νu + grad p = f, div u = 0 in (−1, 1) × (−1, 1) which is discretized by a Chebyshev collocation method. This means that the solution is approximated by global Chebyshev polynomials of degree up to N (cf. [5,22,28]). The resulting saddle point system (1.1) has B = C and A being nonsymmetric and not positive real. In this test, we take five values of N : N = 10, 15, 20, 25 and 30. The corresponding matrices A have orders: 319, 704, 1239, 1924 and 2759, respectively. For each N we use incomplete LU factorizations with five drop tolerances τ : τ = 0.1, 0.01, 0.005, 0.001 and 10−4 for splittings of the (1,1) block A. Analogous to Table 2 in Table 7 we show the outer iteration counts, the total inner iteration counts and the CPU time (in second) for the constraint Schur complement preconditioned GMRES. From Table 7 we can see the same conclusions as those from Table 2.

168

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

Fig. 2 plots the eigenvalues of the constraint Schur complement preconditioned matrix G −1 A arising from N = 15 (the order of G −1 A is 704), τ = 0.1 and τ = 0.01. From Fig. 2 we can see that the eigenvalues corresponding to τ = 0.1 and to τ = 0.01 are both clustered around the point 1. While the eigenvalues corresponding to τ = 0.01 cluster around 1 more tightly than those corresponding to τ = 0.1. Acknowledgements I gratefully acknowledge two reviewers and my editor Dr. Eric de Sturler for many comments and suggestions that helped improve this paper greatly. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

Z.-Z. Bai, M.K. Ng, On inexact preconditioners for nonsymmetric matrices, SIAM J. Sci. Comput. 26 (2005) 1710–1724. R.E. Bank, B.D. Welfert, H. Yserentant, A class of iterative methods for solving saddle point problems, Numer. Math. 56 (1990) 645–666. M. Benzi, G.H. Golub, A preconditioner for generalized saddle point problems, SIAM J. Matrix Anal. Appl. 26 (2004) 20–41. M. Benzi, G.H. Golub, J. Liesen, Numerical solution of saddle point problems, Acta Numer. 14 (2005) 1–137. C. Bernardi, C. Canuto, Y. Maday, Generalized inf-sup conditions for Chebyshev spectral approximation of the Stokes problem, SIAM J. Numer. Anal. 25 (1998) 1237–1271. M.A. Botchev, G.H. Golub, A class of nonsymmetric preconditioners for saddle point problems, SIAM J. Matrix Anal. Appl. 27 (2006) 1125–1140. A. Bouras, V. Frayssé, Inexact matrix-vector products in Krylov methods for solving linear systems, SIAM J. Matrix Anal. Appl. 26 (2005) 660–678. J.H. Bramble, J.E. Pasciak, A.T. Vassilev, Analysis of the inexact Uzawa algorithm for saddle point problems, SIAM J. Numer. Anal. 34 (1997) 1072–1092. J.H. Bramble, J.E. Pasciak, A.T. Vassilev, Uzawa type algorithms for nonsymmetric saddle point problems, Math. Comp. 69 (2000) 667–689. F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. Z.-H. Cao, Fast Uzawa algorithm for generalized saddle point problems, Appl. Numer. Math. 46 (2003) 157–171. Z.-H. Cao, Fast Uzawa algorithms for solving non-symmetric stabilized saddle point problems, Numer. Linear Algebra Appl. 11 (2004) 1–24. Z.-H. Cao, A note on constraint preconditioning for nonsymmetric indefinite matrices, SIAM J. Matrix Anal. Appl. 24 (2002) 121–125. E. de Sturler, J. Liesen, Block-diagonal and constraint preconditioners for nonsymmetric indefinite linear systems. Part I: Theory, SIAM J. Sci. Comput. 26 (2005) 1598–1619. H.S. Dollar, Constraint-style preconditioners for regularized saddle point problems, SIAM J. Matrix Anal. Appl. 29 (2007) 672–684. H.S. Dollar, A.J. Wathen, Approximate factorization constraint preconditioners for saddle-point matrices, SIAM J. Sci. Comput. 27 (2006) 1555–1572. H.C. Elman, Preconditioning for the steady-state Navier–Stokes equations with low viscosity, SIAM J. Sci. Comput. 20 (1999) 1299–1316. H.C. Elman, G.H. Golub, Inexact and preconditioned Uzawa algorithms for saddle point problems, SIAM J. Numer. Anal. 31 (1994) 1645– 1661. H.C. Elman, D.J. Silvester, Fast nonsymmetric iterations and preconditioning for Navier–Stokes equations, SIAM J. Sci. Comput. 17 (1996) 33–46. G.H. Golub, A.J. Wathen, An iteration for indefinite systems and its application to the Navier–Stokes equations, SIAM J. Sci. Comput. 19 (1998) 530–539. M. Gunzburger, Finite Element Methods for Viscous Incompressible Flows, Academic Press, San Diego, 1989. W. Heinrichs, Splitting techniques for the unsteady Stokes equations, SIAM J. Numer. Anal. 35 (1998) 1646–1662. I.C.F. Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput. 23 (2001) 1050–1051. C. Keller, N.I.M. Gould, A.J. Wathen, Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl. 21 (2000) 1300– 1317. L. Lukšan, J. Vlˇcek, Indefinitely preconditioned inexact Newton method for large sparse equality constrained non-linear programming problems, Numer. Linear Algebra Appl. 5 (1998) 219–247. J.A. Meijerink, H.A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a M-matrix, Math. Comput. 13 (1977) 631–644. M.F. Murphy, G.H. Golub, A.J. Wathen, A note on preconditioning for indefinite linear systems, SIAM J. Sci. Comput. 21 (2000) 1969–1972. A. Quarteroni, A. Valli, Numerical Approximation of Partial Different Equations, second, corrected printing, Springer-Verlag, Berlin, 1997. Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Publishing Company, Boston, 1996. Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM J. Sci. Comput. 14 (1993) 461–469. Y. Saad, M.H. Schultz, GMRES:A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869. C. Siefert, E. de Sturler, Preconditioners for generalized saddle-point problems, SIAM J. Numer. Anal. 44 (2006) 1275–1296. D.J. Silvester, N. Kechkar, Optimal low order finite element methods for incompressible flow, Computer Methods Appl. Mech. Engrg. 111 (1994) 357–368.

Z.-H. Cao / Applied Numerical Mathematics 59 (2009) 151–169

169

[34] D.J. Silvester, Stabilized bilinear-constant velocity-pressure finite elements for the conjugate gradient solution of the Stokes problem, Computer Methods Appl. Mech. Engrg. 79 (1990) 71–86. [35] D.J. Silvester, A.J. Wathen, Fast iterative solution of stabilized Stokes systems II: Using general block preconditioners, SIAM J. Numer. Anal. 31 (1994) 1352–1367. [36] V. Simoncini, D.B. Szyld, Theory of inexact Krylov subspace methods and applications to scientific computing, SIAM J. Sci. Comput. 25 (2003) 454–477. [37] L.N. Trefethen, Spectral Methods in Matlab, SIAM, Philadelphia, 2000. [38] H.A. van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University Press, Cambridge, 2003. [39] A.J. Wathen, D.J. Silvester, Fast iterative solution of stabilized Stokes systems I: Using simple diagonal preconditioners, SIAM J. Numer. Anal. 30 (1993) 630–649.