PRECONDITIONING BY AN EXTENDED MATRIX TECHNIQUE FOR CONVECTION-DIFFUSION-REACTION EQUATIONS

. In this paper we consider a preconditioning technique for the ill-conditioned systems arising from discretisations of nonsymmetric elliptic boundary value problems. The rectangular preconditioning matrix is constructed via the transfer operators between successive discretization levels of the initial problem. In this way we get an extended, square, singular, consistent, but mesh independent well-conditioned linear system. Numerical experiments are presented for a 2D convection-diﬀusion-reaction problem.


Let (1.1)
Ax * = b be a linear system of equations, with A an n × n invertible matrix, b ∈ R n and x * its unique solution.If A is ill-conditioned a wide class of preconditioning techniques has been developed in the last 30 years (see [1], [3], [5], [7], [11] and references therein).Usually these methods are of the form (1.2) where the n × n invertible matrices P and Q are constructed in an appropriate way which ensures a much smaller or independent on n condition number for the preconditioned matrix P AQ in (1.2).In this paper we shall consider a more general preconditioning technique which will include (1.2) as a particular case.For this, let m ≥ n and S be a full row rank n × m matrix.Then, the system (1.1) will be transformed as The m × m matrix Â is no more invertible for m > n (because rank( Â) = n), but the system (1.3) is consistent by construction.Moreover, for any solution x ∈ R m of (1.3), by multiplying on the left of (1.3) with S we get SS t AS x = SS t b ⇔ GAS x = Gb ⇔ A(S x) = b ⇔ S x = x * , where G is the n × n symmetric and positive definite matrix (SPD, for short) defined by (1.5) and the superscript t denotes the transpose of a matrix.For any square matrix B we shall denote by σ * (B) the set of all its nonzero eigenvalues and by cond(B) its spectral condition number defined by With respect to the preconditioning procedure (1.3)-(1.4)we are interested whether it exists a constant c > 0, independent on n and m such that A first answer to this question has been given by M. Griebel in [6].In this respect, for symmetric elliptic boundary value problems (b.v.p.), he used for the construction of the preconditioner S in (1.4) the multigrid transfer operators between successive discretization levels of the initial problem.Unfortunately, Griebel's approach essentially uses the symetry of the b.v.p.In this paper, we extend Griebel's multilevel preconditioning procedure to nonsymmetric elliptic b.v.p.The paper is organized as follows: in section 2 we describe Griebel's preconditioner together with the related results in the case of symmetric eliptic b.v.p.In section 3 we introduce the general algebraic framework which gives us the possibility to extend Griebel's ideas to nonsymmetric b.v.p. and prove that the preconditioned matrix has a mesh independent condition number (Theorem 1).In section 4 we present numerical experiments with two versions of our preconditioning technique for a class of 2D convection-diffusion-reaction problems.

GRIEBEL'S MULTILEVEL PRECONDITIONING
In [5] M. Griebel considered the following boundary value problem where Ω = (0, 1) d , (usually d = 1, 2, 3) and L is a second order elliptic differential operator on Ω.The problem (2.1) can be discretized by appropriate finite differences (see [3]) or finite element techniques (see [7], [9]).In the later case we may also suppose that a variational formulation is available in the following form: find u ∈ U (Ω) with where U (Ω), V (Ω) are Hilbert spaces of real valued functions defined on Ω, a : are the L 2 (Ω) scalar product and norm.We shall also suppose that the bilinear functional a is such that the variational formulation (2.2) has a unique solution u ∈ U (Ω).In [5] M. Griebel considered a symmetric elliptic b.v.p. of the form (2.1)-(2.2) for which U (Ω) = V (Ω) = H 1 0 (Ω) and a bilinear, symmetric, bounded and coercive, i.e.
where v 2 Remark 2.1.For g = 0 there is a straightforward way to obtain a variational formulation as (2.2)-(2.3),whereas in the nonhomogeneous case, g = 0, we may use the approach from Chapter 7 in [7].
n k } a standard finite element basis (e.g.piecewise d-linear, see [9]).Then the linear system associated to (2.2) is From the properties of the bilinear functional a it results that the matrix A k is SPD.Let now be a sequence of spaces of piecewise d-linear functions associated to a sequence of uniform, equidistant, nested grids nq } the finite element basis in V q .Let also Bk ⊂ V k and m k be given by The functions from Bk are linearly dependent and generate the subspace We consider the n q+1 × n q grid transfer matrix I q+1 q given by (2.10)I q+1 q ij = c ij and for q = 1, 2, . . ., k − 1 define the n k × n q matrices S k q by (2.11) , in which the last n k × n k block is the unit matrix.We then consider the It results that the m k × m k matrix Âk is symmetric and positive semidefinite.M. Griebel proved (see [5] and [6] references therein) that the spectral condition number of the preconditioned matrix Âk is mesh independent.He's proof is essentially based on the symmetry of A k (and Âk ) and the fact that the n k × n k matrix G k defined as in (1.5), i.e. (2.15) is spectrally equivalent (see the definition (3.1) in section 3) with the inverse of the standard discretized Laplacian ∆ k (5-point stencil, see [7]), which for picewise finite element basis functions and a variational formulation on In the next section we shall present an extension of this procedure for nonsymmetric problems as (2.1).

THE GENERAL ALGEBRAIC FORMULATION
We shall start the presentation of our extension by a brief replay of some results related to spectrally equivalent matrices.In this sense we introduce the following notations: for an n × n SPD matrix A we shall denote by σ(A), ρ(A) its spectrum and spectral radius and by λ min (A), λ max (A) its smallest and biggest eigenvalue.For an arbitrary m × n matrix T , N (T ), R(T ) will denote its null space and range, respectively; the notations •, • , • will be used for the euclidean scalar product and norm on some vector space R q .Definition 3.1.For A and B two n × n SPD matrices, we shall say that A is spectrally equivalent with B if there exist two positive constants α 1 , α 2 > 0, independent on n such that If we write this by A ≈ B the following results are known (see e.g.[1,3,4]).
(i) The matrix A is spectrally equivalent with B if and only if for a decomposition of the form B = CC t we have for some a 1 , a 2 independent on n.Moreover, inequalities (3.2) are independent on the decomposition of B.
2) makes the connection between the spectral equivalence and preconditioning.Indeed, by taking into account that the spectral condition number of an SPD matrix T is defined as (see also (1.6)) (3.3) cond(T ) = λmax(T ) λ min (T ) , from (3.2) it results that, if the matrices A and B are spectrally equivalent, then B is a "good preconditioner" for A, i.e. (3.4) cond(C −1 AC −t ) ≤ α 2 α 1 .let us know consider a general (nonsymmetric) system of the form (1.1) and let M, R be the matrices defined by it results that both matrices are normal (see e.g.[1]).We shall suppose that M is positive definite, i.e.
Our extension of the previous Griebel's preconditioning procedure is based on the following assumptions.
The following results shows that the preconditioned matrix Â from (1.4) is well-conditioned.
Theorem 3.4.In the above hypothesis and under the assumptions A1 and A2 we have for the matrix Â in For proving Theorem 1 we need two auxiliary results which will be presented in what follows.
Lemma 3.5.We have the equality Proof.Let first λ ∈ σ * ( Â).Then, for some nonzero vector z we have the following sequence of equalities, in which the last one shows that λ ∈ σ(GA).
Conversely, let λ ∈ σ(GA).Then, λ = 0 and for some nonzero vector y we have GAy = λy.Using the invertibility of the SPD matrix G we define w = G −1 y = 0 and get GAGw = λGw ⇐⇒ SS t ASS t w = λSS t w.
Thus, because the application S t : R n −→ R m is injective, the vector z = S t w will be nonzero and S Âz = λSz, i.e. (3.12) Âz − λz ∈ N (S).
Lemma 3.6.Let G = CC t with C an n × n invertible matrix and Ā, M , R defined by with A, M, R from (3.5).Then Proof.Because M (see ( ) is SPD, so will be M from (3.13).Then, by also using Assumption A1 we get λ min ( M ) = min Moreover, from (3.13) and (3.9) it follows that Now, by using (3.11) and some well-known properties of the spectrum and spectral radius of matrices, we successively get Then, from (3.15), (3.17), (1.6), (3.14) and (3.9), we obtain (3.10) and the proof is complete.
Remark 3.7.We have to observe that, beside the specific properties mentioned in the above assumptions A1 and A2, our preconditioning method requests only the invertibility of the system matrix in (1.1) and the fact that its symmetric part is SPD.These properties are fulfilled by a large class of finite element or finite differences discretizations of elliptic boundary value problems (b.v.p.), see e.g.[8]).

NUMERICAL EXPERIMENTS
We considered in our experiments the convection-diffusion-reaction problem (4.1) where γ ∈ (0, ∞) and the right hand side f such that the unique exact solution is u(x, y) = xy(1 − x)(1 − y)e xy (see [3]).The problem was discretized using the classical 5-point stencil finite differences on k successive grids.On each grid, the number of interior nodes is n q = (2 q − 1) 2 , q = k, . . ., 1, such that the system on the finest level (2.4) has dimension n k .We used in our experiments two types of intergrid transfer operators, I q+1 q and J q+1 q , q = k − 1, . . ., 1 defined (in stencil notation) by (4.2) The preconditioning n k × m k matrix S k was constructed as in (2.11)-(2.12)and the extended preconditioned system and matrix G k as in (2.13)-(2.15).For solving the initial system (2.4) and the preconditioned one (2.13)-(2.14)we used the CGLS (Conjugate Gradient for Least Squares) algorithm from [1], with the stopping rules, respectively The results (number of iterations versus the dimension of the finest grid n k and the convection coefficient γ (for the reaction coefficient δ = 0)) are presented in Tables 1-3.We observe in Tables 2-3 the mesh independence behaviour of the preconditioned system for both choices of intergrid transfer operators, together with γ-dependence for fixed k and n k .Moreover in Table 4 we can see that Âk remains a sparse matrix and (because of the mesh independence of the preconditioning) the computational time for solving the system (2.13)-(2.14)(constructed with I q+1 q from (4.2) is much less, for bigger dimensions, than for solving the nonpreconditioned one (2.4).

CONCLUSIONS AND FURTHER WORK
In this paper we presented an extension of Griebel's preconditioning technique from [5], for symmetric elliptic boundary value problems to nonsymmetric problems if the form (4.1).The extension, formulated in Section 3 in a very general algebraic form is based on the assumptions A1 and A2, which controls the symmetric and, respectively skew-symmetric part of the system matrix A. We applied this extension to the (nonsymmetric) boundary value problem (4.1), using a finite differences discretization and a classical coarsening.We used the two versions of intergrid transfer operators from (4.2).In both cases, the results presented in Tables 2 and 3 show a mesh-independence behaviour.For I k k−1 the symmetric part M k of A k is exactly the 5-point stencil Laplacian ∆ k (see [3], [7]), thus the assumption A1 holds.Unfortunately we have not yet a similar result for the J k k−1 operators and also for the second assumption A2, in both cases from (4.2).Work is in progress on this direction.
t AS, b = S t b.

Table 4 .
Sparsity and computational time (for J q+1 k spa(A) spa( A k ) T ime(A k ; b k ) T ime( A k ; b