Caterina Calgaro
email: caterina.calgaro@univ-lille1.fr

Jean-Paul Chehab
email: jean-paul.chehab@u-picardie.fr

Yousef Saad
email: saad@cs.umn.edu

Incremental Incomplete LU factorizations with applications

Keywords: Preconditioning, Incomplete LU factorization, incremental LU

This paper addresses the problem of computing preconditioners for solving linear systems of equations with a sequence of slowly varying matrices. This problem arises in many important applications. For example, a common situation in computational fluid dynamics, is when the equations change only slightly, possibly in some parts of the physical domain. In such situations it is wasteful to recompute entirely any LU or ILU factorizations computed for the previous coefficient matrix. A number of techniques for computing incremental ILU factorizations are examined. For example we consider methods based on approximate inverses as well as alternating techniques for updating the factors L and U of the factorization.

Introduction

When solving large sparse linear systems of equations, it is generally accepted that preconditioned Krylov subspace methods are effective alternatives to direct solution methods which offer a good compromise between efficiency and robustness. These methods combine an accelerator such as GMRES or Bi-CGSTAB and a preconditioner such the ILU factorization with threshold or with level-of-fill, see, e.g., [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Axelsson | Iterative Solution Methods[END_REF] for details. Here, a few techniques are presented for constructing ILU-type preconditioners in an incremental manner, i.e., starting from an already computed approximate factorization. There are many reasons why such techniques may be desirable.

For example, a common problem which arises in many complex applications is to solve a sequence of linear systems of the form:

A k x k = b k , (1)
for k = 1, 2, ... In these applications A k does not generally change too much from one step k to the next, as it is often the result of a continuous process (e.g. A k can represent the discretization of some problem at time t k ). The problem then is to solve each consecutive system effectively, by taking advantage of earlier systems if possible. This problem arises for example in computational fluid dynamics, when the equations change only slightly possibly in some parts of the domain. A similar situation was recently considered in [START_REF] Birken | Preconditioner updates applied to CFD model problems[END_REF].

There are many possible other applications of incremental ILU factorizations. One that is often neglected concerns the construction of preconditioners progressively within the whole iterative cycle. Thus, when solving a linear system we do not know in advance how accurate the preconditioner must be for the iteration to converge fast enough. One could then start with a rough factorization which is gradually improved within a Krylov subspace acceleration. As soon as there are signs of stagnation -or slow convergence -one can attempt to improve the current factors by running one step of improvement with the help of one of several procedures that will be introduced in this paper.

It is also possible to devise strategies to build ILU factorizations in this way. Starting with the SSOR(1) factorization, one can compute successive improvements until a satisfactory factorization is obtained. As will be seen, a by-product of one of the techniques introduced is that one can compute an exact LU factorization by a sort if alternating procedure. We will establish the surprising result that this iterative procedure in fact converges in N steps, where N is the matrix size.

A somewhat related problem which was examined in the past is that of updating the (exact) LU factorization of a given matrix when a very small number of its entries change. This particular case received much attention in the power systems literature, see, e.g., [START_REF] Chan | Partial matrix refactorization[END_REF]. In contrast, the problem of updating an existing ILU factorization has received little attention so far. Birken et al. [START_REF] Birken | Preconditioner updates applied to CFD model problems[END_REF], consider this problem specifically for CFD equations and proposed an inexpensive approach to update one of the factors whey they undergo small changes. The basic idea of their method is quite simple and is discussed in Section 4. [START_REF] Benzi | A sparse approximate inverse preconditioner for the conjugate gradient method[END_REF].

Updating preconditioners is quite easy when approximate inverse preconditioners are used because these methods typically involve explicit updates. We will start by exploring techniques for improving incomplete LU factorizations by exploiting approximate inverse strategies. We will then discuss a set of techniques based on alternating procedures for correcting an ILU factorization.

In the remainder of this section we state the specific problem addressed in this paper, and outline in broad terms the two classes of methods considered for solving it. Some notation which will be used throughout the paper will also be introduced.

The problem

Suppose we have an approximate factorization which we write in the form

A = LU + R, ( 2
)
where R is the error in the factorization, which corresponds to the matrix of the terms that are dropped during the factorization procedure. Our primary goal is to improve this factorization, i.e., to find sparse matrices X L , X U such that L + X L , U + X U is a better pair of factors than L, U . Consider the new error A -LU for the new factors.

A

-(L + X L )(U + X U ) = (A -LU ) -X L U -LX U -X L X U . (3)
Ideally, we would like to make the right-hand side equal to zero. Denoting by R the matrix A -LU , this means we would like to solve:

X L U + LX U + X L X U -R = 0 , (4)
for X L and X U . The above equation, which can be viewed as a form of Riccati equation, is nonlinear in the pair of unknowns X L , X U . In the symmetric case, where LU is replaced by the Cholesky factorization, and where we set X U = X T L ≡ X by symmetry, then (4) becomes

XL T + LX T + XX T -R = 0 , (5)
the standard Algebraic Riccati Equation.

In our case, we only consider (4) which we would like to solve approximately. For this we begin by neglecting the quadratic term X L X U , leading to:

X L U + LX U -R = 0 . (6)
We will consider two broad paths to tackle the above problem. First, we will consider a few approaches in the spirit of approximate inverse-type techniques which will try to find sparse triangular factors X L , X U that approximately minimize the Frobenius norm of the left-hand side of [START_REF] Cahouet | Some fast 3D finite element solvers for the generalized Stokes problem[END_REF]. Second, we exploit an alternating procedure at the matrix level which fixes U (i.e., we set X U = 0) and solves for X L , and then fixes the resulting L and solve for X U , and repeat the process. It is also possible to use a similar alternating procedure but at the vector levelwhere columns of L and rows of U are updated alternatively. This approach, which is similar to a well-known one used for solving Lyapunov equations, is not considered here. Similarly, methods based on differential equations similar to those discussed in [START_REF] Chehab | Differential equations and solution of linear systems[END_REF][START_REF] Chehab | Matrix differential equations and inverse preconditioners[END_REF] can also be of interest in this context but will not be considered in this paper.

Notation

We will often need to extract lower and upper triangular parts of matrices. Given an N × N matrix X, we denote by X ⌊ \ the strict lower triangular part of X, and by X \ ⌉ the upper triangular part of X (which includes the diagonal of X). The unit lower triangular matrix obtained by taking the lower triangular part of X and replacing its diagonal entries by ones, is denoted by X ⌊ . In other words X ⌊ = I + X ⌊ \ , the identity matrix plus the strict lower part of X. Note the following easy to verify properties of these operators

(X + Y ) ⌊ \ = X ⌊ \ + Y ⌊ \ ; (X + Y ) \ ⌉ = (X + Y ) \ ⌉ ; (7) (X ⌊ Y ⌊ \ ) \ ⌉ = 0; (X ⌊ \ Y ⌊ ) \ ⌉ = 0; (8) (X \ ⌉ Y \ ⌉ ) ⌊ \ = 0 . (9)
The fact that a matrix is upper triangular can be conveniently expressed by the equality X \ ⌉ = X. Similarly, a matrix is strict lower triangular when X ⌊ \ = X and is unit lower triangular when X ⌊ = X. This notation is specifically designed for the LU factorization which involves unit lower triangular and upper triangular matrices. The inner product of two N × N matrices, is defined as:

X, Y = Tr(Y T X). (10)
The Frobenius norm . F is the 2-norm associated with this inner product, i.e., X F = [Tr(X T X)] 1/2 .

Techniques based on approximate inverses

Standard approximate inverse techniques attempt to find an approximate inverse to A which is sparse, see, e.g., [START_REF] Benzi | A sparse approximate inverse preconditioner for the conjugate gradient method[END_REF][START_REF] Benzi | Stabilized and block approximate inverse preconditioners for problems in solid and structural mechanics[END_REF][START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear systems[END_REF][START_REF] Kolotinina | Incomplete block factorizations as preconditioners for sparse spd matrices, in Recent advances in iterative methods[END_REF][START_REF] Kolotilina | On a family of two-level preconditionings of the incomplete block factorization type[END_REF][START_REF]Factorized sparse approximate inverse preconditionings I.Theory[END_REF][START_REF] Grote | Parallel preconditioning and approximate inverses on the connection machine[END_REF][START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF]. One such approach used in APINV (see, e.g., [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF]) is to approximately minimize the Frobenius norm,

I -AM F (11)
over sparse matrices M which have a pre-specified pattern. We are interested in a by-product of this general class of methods which we term "Sparse Matrix Corrections".

Sparse Matrix Corrections

If B is a given preconditioner to A, we can seek a matrix M such that AM approximates B instead of the identity matrix as is done in [START_REF] Chehab | Differential equations and solution of linear systems[END_REF]. The approximation can be sought in the Frobenius norm sense. More precisely, let B = LU be a given ILU preconditioner of A. We can seek a matrix M such that AM approximates B = LU by trying to approximately minimize B -AM F . Note that B = LU is generally a sparse matrix. The corrected preconditioning operation would then become M B -1 = M U -1 L -1 . This approach was first advocated in [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF] (see also [28, sec. 10.5.8]), and was also exploited in [START_REF] Chow | Approximate inverse techniques for block-partitioned matrices[END_REF]. It was later discussed in [START_REF] Holland | Sparse approximate inverses and target matrices[END_REF] where the term "target matrices" was introduced for the matrix B. In the above equation M can be obtained by approximately solving the linear systems Am j = b j , where b j , the j-th column of B is sparse, in order to find the j-th (sparse) column m j of M .

An alternative to solving individual successive linear systems Am j = e j , is to seek a global iteration, i.e., an iteration which updates the whole matrix at once as in

M new = M + αG (12)
where G, the 'search direction', is a certain matrix. It is possible to apply a form of the steepest descent algorithm to try to minimize

F (M ) = B -AM 2 F , (13)
which is a quadratic function on the space of N × N matrices, viewed as objects in R N 2 . This requires the gradient of the objective function of F (M ) in [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF]. The natural inner product on the space of matrices, with which the function ( 13) is associated, is the inner product [START_REF] Chehab | Matrix differential equations and inverse preconditioners[END_REF].

In descent-type algorithms the scalar α in ( 12) is selected to minimize the objective function associated with M new . The simplest choice for the descent direction G is to take it to be the residual matrix R = B-AM , where M is the new iterate. The corresponding descent algorithm is called the "Minimal Residual (MR)" algorithm. It is easy to see that R is indeed a descent direction for F . See [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF] for details on the resulting global Minimal Residual algorithm.

Another popular choice is to take G to be the direction of steepest descent, i.e., the direction opposite to the gradient of F . A result of [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF] can be expanded to show that the array representation of the gradient of F with respect to M is the matrix G = -2A T R. As noted in [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF] this option can lead to a very slow convergence in some cases because it is essentially a steepest descent-type algorithm applied to the normal equations. In either global steepest descent or minimal residual, we need to form and store the G matrix explicitly.

The tools just discussed can be adapted to develop algorithms for improving a given ILU factorization, by alternatively correcting the L and U factors in turn.

Alternating L-U descent methods

Consider first equation ( 4) with the particular choice X L = 0, which corresponds to updating U while L is kept frozen. Then, a sparse X U is sought which minimizes,

F (X U ) = A -L(U + X U ) 2 F = R -LX U 2 F ( 14
)
where we set R ≡ A -LU , the current factorization error. Similar to the situation of approximate inverses, the optimum X U is trivially known, in this case it equals L -1 R, and this will be exploited in the next section. Here, we seek an approximation only to this exact solution. Note that ( 14) is of the same form as [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF]. We have

R -LX U 2 F = Tr ([R -LX U ] T [R -LX U ]) = R 2 F -2Tr(R T LX U ) + LX U 2 F . Therefore, the gradient of F (X U ) at X U = 0 written in matrix form is simply G = -2 L T R.
What this means is that when X U is a small multiple of L T R then F (X U ) will decrease. With this, a steepest descent method can now be devised. The only issue is that G = L T R is not necessarily upper triangular, so X U which is a multiple of G in this situation will not be upper triangular as desired. We can define G to be the upper triangular part of L T R, i.e., using the notation defined in Section 1.2,

G = [L T R] \ ⌉ , (15)
which represents the upper triangular part of L T R (including the diagonal). We may need to sparsify G matrix even further to limit fill-in. Remarkably, the resulting matrix remains a descent direction under these changes. Here descent direction is to be interpreted as to mean 'non-ascent'.

Lemma 1. The matrix G = [L T R] \ ⌉ and any sparsified version G of it is a descent direction for F at X U = 0.
Proof. Let us denote by X the matrix L T R. Since the gradient is -2X, a matrix G is a descent direction if < -2X, G >≤ 0, or < X, G >= Tr(G T X) ≥ 0. Let η ij be the entries of X and consider any diagonal entry of G T X, where G is the upper triangular part of X. We have:

[G T X] jj = j i=1 η 2 ij ≥ 0 .
Therefore the trace is ≥ 0. If G is further sparsified to G, then the above sum is replaced by a sum over the sparsity pattern of the j-th row of G T , i.e., the sparsified pattern of the j-th column of G, and it remains non-negative.

To determine the best α in [START_REF] Chow | Approximate inverse techniques for block-partitioned matrices[END_REF] we must minimize the quadratic form:

R -αLG 2 F = Tr ([R -αLG] T [R -αLG]) = R 2 F -2αTr(R T LG) + α 2 LG 2 F , yielding α U = R, LG LG, LG = Tr(R T LG) Tr ((LG) T LG) . ( 16
)
Formulas for the L part can be readily obtained by repeating the above argument, replacing the objective function [START_REF] Davis | Direct methods for sparse linear systems[END_REF] by

A -(L + X L )U 2 F = R -X L U 2 F .
The array representation of the gradient becomes -2RU T . We seek a correction X L to L that is strict lower triangular, so we define G L = [RU T ] ⌊ \ and

α L = R, G L U G L U, G L U = Tr(R T G L U ) Tr ((G L U ) T G L U ) . ( 17
)
Note that the denominators in ( 16) and ( 17) may be computed as LG 2 F and G L U 2

F

respectively. An alternating algorithm which updates U and L in turn, can now be sketched.

Algorithm 1 (MERLU, Minimal Energy Residual descent for LU )

1.

Select an initial pair L, U (with

L = L ⌊ , U = U \ ⌉ ) 2.
Until convergence Do 3.

Compute

R := A -LU 4. Compute G = [L T R] \ ⌉ 5.
Apply numerical dropping to G 6.

Compute α = R, C / C 2 F , where C = LG.

7.

Compute U := U + αG 8.

Compute R := A -LU 9.

Compute G := [RU T ] ⌊ \ 10.

Apply numerical dropping to G 11.

Compute α = R, C / C 2 F where C = GU . 12.

Compute L := L + αG 13. EndDo With the above approach it is guaranteed that the residual norm R F will not increase. A convenient feature of the method is that it can tolerate dropping as indicated by Lemma 1. One possible practical concern is the cost of this procedure. If we perform one iteration to improve a given preconditioner, the computation will be dominated by the formation of the gradient matrices G in Lines 4 and 9 and the calculation of the scalars α in lines 6 and 11. These involve the products of two sparse matrices. Note that the matrix C need not be stored. It can be calculated one column at a time, for example. Then the norm of the column and the inner product of this column with the corresponding column of R can be evaluated to compute parts of R, C and C, C and then the column can be discarded.

Alternating Sparse-Sparse iterations

Consider equation (4) in which we set X U = 0. If U is nonsingular we obtain:

X L U = R → X L = RU -1 . (18)
Thus, the correction to L can be obtained by solving a sparse triangular linear system with a sparse right-hand side matrix, i.e., the system U T X T L = R T . However, as was noted before, the updated matrix L + X L obtained in this way is not necessarily unit lower triangular and it is therefore required to apply the [ . ] ⌊ operation for L + X L or the [ . ] ⌊ \ operation for X L . This procedure can be repeated by freezing U and updating L and vice-versa alternatingly, leading to an iteration which can be summarized by the following two equations:

U k+1 = U k + [L -1 k (A -L k U k )] \ ⌉ , (19)
L k+1 = L k + [(A -L k U k+1 )U -1 k+1 ] ⌊ \ . (20)
A feature of equations [START_REF] Golub | Matrix Computations[END_REF][START_REF] Grote | Parallel preconditioning and approximate inverses on the connection machine[END_REF] is that they simplify. For [START_REF] Golub | Matrix Computations[END_REF] for example, we have

U k+1 = U k + [L -1 k (A -L k U k )] \ ⌉ = [U k + L -1 k (A -L k U k )] \ ⌉ = [L -1 k A] \ ⌉ .
A similar expression for L k+1 is easily derived and this leads to an iteration of the form:

U k+1 = [L -1 k A] \ ⌉ , (21)
L k+1 = [AU -1 k+1 ] ⌊ . (22)
Note that the above two formulas are not intended for practical use. It is generally more appealing to use [START_REF] Golub | Matrix Computations[END_REF][START_REF] Grote | Parallel preconditioning and approximate inverses on the connection machine[END_REF] because one can sparsify A -L k U k (say) keeping only the largest entries and making the solve much less inexpensive. One issue in the above equations is the fact that the algorithm will break down when U k+1 is singular, and this may indeed happen, as the following example shows.

Example

A =   1 1 1 1 2 2 1 2 4   → L 0 =   1 0 0 2 1 0 -3 -1 1   → L -1 0 A =   1 1 1 -1 0 0 3 5 7   → U 1 =   1 1 1 0 0 0 0 0 7   .
This issue is not too different from that encountered in incomplete LU factorizations with threshold. In fact any form of LU factorization implicitly involves solves with lower and upper triangular matrices and when dropping is applied, these matrices can become singular. The occurrence of this breakdown is rare in practice. Note that at the limit, U k will approach the U factor of the LU factorization of A (as will be shown in Proposition 1), a nonsingular matrix.

The resulting alternating algorithm is given below.

Algorithm 2 (ITALU, Iterative Threshold Alternating Lower-Upper correction)

1. Given: A, U 0 , L 0 (with U 0 = [U 0 ] \ ⌉ ; L 0 = [L 0 ] ⌊ ) 2. For k = 0, • • • , Do: 3. Compute R k = A -L k U k 4. Compute X U = [L -1 k R k ] \ ⌉ 5.
Apply numerical dropping to X U 6.

U k+1 = U k + X U 7. If det(U k+1 ) == 0 Abort "Singular U reached" 8. Compute R k+1/2 = A -L k U k+1 9. Compute X L = [R k+1/2 U -1 k+1 ] ⌊ \ 10.
Apply numerical dropping to X L 11.

L k+1 = L k + X L 12. EndDo
It may be thought that the triangular systems in Lines 4 and 9 of the algorithm are expensive to solve. However, one must remember that the right-hand side is sparse and there exist effective techniques for this task. This is discussed in detail in Section 4. A variation to the dropping steps in Lines 5 and 10, would be to drop terms from the residual matrices R k and R k+1/2 instead of the matrices X U and X L .

Expressions [START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF] can be further simplified by replacing A by one of its triangular parts. This is a result of the following simple lemma which will also be helpful in the next section.

Lemma 2. Let X be any square matrix and let P be a unit lower triangular matrix and Q an upper triangular matrix. Then the following equalities hold:

[P X] \ ⌉ = [P X \ ⌉ ] \ ⌉ ; [XQ] ⌊ \ = [X ⌊ \ Q] ⌊ \ . ( 23
)
Proof. The first relation results from the equalities ( 7) and ( 8). Indeed,

[P X] \ ⌉ = [P (X ⌊ \ + X \ ⌉ )] \ ⌉ = [P X ⌊ \ ] \ ⌉ + [P X \ ⌉ )] \ ⌉ = 0 + [P X \ ⌉ )] \ ⌉ .
The second equality can be shown similarly from ( 7) and ( 9):

[XQ] ⌊ \ = [(X ⌊ \ + X \ ⌉ )Q] ⌊ \ = [X ⌊ \ Q] ⌊ \ + [X \ ⌉ Q] ⌊ \ = [X ⌊ \ Q] ⌊ \ + 0.
Note that other similar equalities can also be shown. For example we will later exploit the relation:

[XQ] ⌊ \ = [X ⌊ Q] ⌊ \ .
Using the lemma with

P = L -1 k , Q = U -1 k+1
, and X = A, equations [START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF] become:

U k+1 = [L -1 k A \ ⌉ ] \ ⌉ , (24)
L k+1 = [A ⌊ \ U -1 k+1 ] ⌊ . (25)

Convergence results in the dense filter case

Filtering consists of keeping only specific wanted entries in a matrix, and dropping the others, i.e., replacing them with zeros. For example, applying an upper triangular filter to a matrix M amounts to only keeping its upper triangular part. An equivalent term often used in the literature is that of a 'mask'. The filter matrix can be dense or sparse. Next, we establish a few properties of the sequences U k and L k in the particular case where we use dense triangular filters. This only means that the iterates L k , U k are evaluated by the expressions [START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF] exactly (no other dropping is performed). We start with the following lemma.

Lemma 3. Let two square matrices P and Q be unit lower triangular and upper triangular respectively, with Q nonsingular. Then for any square matrix M we have:

[M ] ⌊ = [[M Q -1 ] ⌊ Q] ⌊ and [M ] \ ⌉ = [P [P -1 M ] \ ⌉ ] \ ⌉ . (26)
Proof. The second part of the equation in ( 23) applied to

X = M Q -1 together with the definition [.] ⌊ = I + [.] ⌊ \ , yields immediately [[M Q -1 ] ⌊ Q] ⌊ = [(M Q -1 )Q] ⌊ = [M ] ⌊ .
The second relation in [START_REF] Kolotinina | Incomplete block factorizations as preconditioners for sparse spd matrices, in Recent advances in iterative methods[END_REF] follows similarly by taking X = P -1 M in the first equality of equation ( 23)

We can now prove the following proposition

Copyright c 2009 John Wiley & Sons, Ltd.

Numer. Linear Algebra Appl. 2009; 0:1-1 Prepared using nlaauth.cls Proposition 1. Let A be a matrix that admits the LU decomposition A = LU , where L is unit lower triangular and U is upper triangular. Then, the sequences U k and L k defined by [START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF] and [START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF] starting with an arbitrary (unit lower triangular) L 0 , satisfy the following properties

(i) [A -L k U k+1 ] \ ⌉ = 0, [A -L k+1 U k+1 ] ⌊ \ = 0, k ≥ 0. (ii) If the sequences L k and U k converge, then lim k→+∞ L k = L and lim k→+∞ U k = U . (iii) The algorithm converges in one step if L 0 = L.
Proof. We first prove assertion (i). From [START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF]

it follows that L k U k+1 = L k [L -1 k A] \ ⌉ .
Applying an upper triangular filter to both sides, we obtain, thanks to the second part of Equation ( 26)

[L k U k+1 ] \ ⌉ = [A] \ ⌉ .
This proves the first half of (i). To prove the second half first observe that from [START_REF] Holland | Sparse approximate inverses and target matrices[END_REF] we have

[L k+1 U k+1 ] ⌊ \ = [(L k+1 ) ⌊ \ U k+1 ] ⌊ \ . Then, using the relation L k+1 = [AU -1
k+1 ] ⌊ , and invoking Equation ( 26) once more yields:

[L k+1 U k+1 ] ⌊ \ = [[AU -1 k+1 ] ⌊ U k+1 ] ⌊ \ = [A] ⌊ \ . To prove assertion (ii), we introduce L = lim k→+∞ L k and Ū = lim k→+∞ U k . Taking limits in (i) yields the relations [ L Ū ] \ ⌉ = [A] \ ⌉ and [ L Ū ] ⌊ \ = [A] ⌊ \ .
Hence L Ū = A and the uniqueness of the LU factorization implies that L = L and Ū = U . Finally, when we take

L 0 = L, then U 1 = [L -1 LU ] \ ⌉ = [U ] \ ⌉ = U , proving (iii).
In fact, if one of the sequences L k or U k converges, then both of them will converge to their limits L and U respectively. This follows from the following relations

U k+1 -U = [L -1 k A -U ] \ ⌉ = [L -1 k A -L -1 A] \ ⌉ = [(L -1 k -L -1 )A] \ ⌉ L k+1 -L = [AU -1 k+1 -L] ⌊ \ = [AU -1 k+1 -AU -1 ] ⌊ \ = [A(U -1 k+1 -U -1
)] ⌊ \ . These two relations show that as long as there is no break-down (i.e., U k+1 is nonsingular each time it is computed from ( 24)) the sequences L k , U k are well defined, and they converge hand-in hand, in the sense that when L k converges, then U k also converges and vice versa.

We can now prove the following convergence result which establishes the surprising fact that under the condition that the filters are triangular and dense, the sequence of L k 's and U k 's terminates in a finite number of steps.

Theorem 1. Let an initial L 0 be given such that the sequences L k , U k are defined for k = 1, 2, • • • , N . Then the sequences U k and L k converge to U and L respectively, in no more than N steps.

Proof. We proceed by induction on the dimension N of the matrices. The result is trivial to establish for the case

N = 1. Indeed, if A = [a], the trivial LU factorization is L = [1], U = [a]
and the sequence L k , U k reaches this factorization in one step.

We now assume that the algorithm converges in no more than N -1 steps for matrices of dimension N -1. We use the notation We also consider the block LU factorization of A,

A = Ã b c d , L k = Lk 0 α k 1 , U k = Ũk γ k 0 δ k . Copyright c
A = L 0 c Ũ -1 1 L Ũ L-1 b 0 δ U with δ = d -c Ũ -1 L-1 b. ( 27
)
The recurrence relations [START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF] of the algorithm give,

U k+1 = Ũk+1 γ k+1 0 δ k+1 = [ L-1 k Ã] \ ⌉ L-1 k b 0 -α k L-1 k b + d . ( 28
)
and

L k+1 = Lk+1 0 α k+1 1 = [ Ã Ũ -1 k+1 ] ⌊ 0 c Ũ -1 k+1 1 . ( 29
)
This results in the relations

Ũk+1 = [ L-1 k Ã] \ ⌉ , Lk+1 = [ Ã Ũ -1 k+1 ] ⌊ , (30)
and

α k+1 = c Ũ -1 k+1 , γ k+1 = L-1 k b, δ k+1 = d -α k L-1 k b. (31)
Now by assumption, the relations [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] show that the sequence of the pairs Lk , Ũk converges in no more than N -1 steps to the L, U factors of Ã, i.e., LN-1 = L and ŨN- [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] and comparing the result with the U part of the factorization [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] we obtain immediately U N = U . Proceeding the same way by taking k = N -1 in [START_REF] Silvester | Efficient preconditioning of the linearized Navier-Stokes equations for incompressible flow[END_REF] and comparing the result with the L factor in [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] shows that L N = L and completes the proof.

1 = Ũ . It follows that LN = L, ŨN = Ũ , α N = c Ũ -1 , γ N = L-1 b, δ N = d -c Ũ -1 L-1 b. Taking k = N -1 in
It is interesting to observe that the above theorem provides an alternative way of computing the (exact) LU factorization of a matrix. Like Gaussian elimination, it is a direct method, though a costlier one in the dense case. Specifically, the full N steps of the algorithm require O(N 4 ) operations to complete in the dense case. However, in general the algorithm is likely to converge in fewer than N steps.

Convergence in the sparse filter case

Consider now the situation where a general sparse filter is applied. This means that a sparsity pattern F k is selected at each step k and the matrices L k , U k of the factorization are restricted to having the pattern of the lower and upper parts of F k , respectively, by dropping any entries outside of this pattern. Recall the notation with filters. The filter F can be considered as an N × N matrix of zeros and ones. A zero entry corresponds to an unwanted term in the pattern, so a nonzero entry in this location is dropped. We denote by X ⊙F the Hadamard product (i.e., component-wise product), of the matrices X and F. The algorithm with a general constant sparse filter would replace equations [START_REF] Grote | Parallel preconditionings with sparse approximate inverses[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF] by

U k+1 = [(L -1 k A) ⊙ F] \ ⌉ = [(L -1 k A)] \ ⌉ ⊙ F (32) L k+1 = [(AU -1 k+1 ) ⊙ F] ⌊ = [(AU -1 k+1 )] ⌊ ⊙ F . ( 33
)
Copyright c 2009 John Wiley & Sons, Ltd.

Numer. Linear Algebra Appl. 2009; 0:1-1 Prepared using nlaauth.cls

The algorithm breaks down immediately if the filter has a zero diagonal entry. These singular filters and excluded from consideration, i.e., all filters considered are such that F ii = 1, ∀i.

A desirable extension of the convergence result of the previous section would be to show that the sequences L k , U k converge to the sparse LU factors L and U of A respectively. For this to be possible, the patterns of the filter matrices F k must contain those of the L and U parts respectively at each step. In fact under the assumption that F k contains the symbolic pattern of Gaussian elimination, the same result as that of the previous theorem holds. Here, by symbolic pattern, we mean the pattern of Gaussian elimination obtained from the graph, i.e., in the situation where exact cancellations do not take place in Gaussian elimination [START_REF] Davis | Direct methods for sparse linear systems[END_REF]. The proof of this result is an immediate extension of that of the previous theorem. The result is stated without proof as a corollary.

Corollary 1. Assume that each of the filter matrices F k , k ≥ 0 contains the (symbolic) sparsity patterns of L and U respectively and let an initial L 0 be given such that L 0 ⊙ F 0 = L 0 and the sequences L k and U k exist for k = 1, • • • , N . Then the sequences U k and L k converge to U and L, respectively, in no more than N steps.

It is interesting to consider how the results of Proposition 1 change when a general sparse filter is applied. Consider the iteration (32), (33), in which the filter F is now constant. We can rewrite the iteration with the help of the upper and lower parts of the filter as follows:

U k+1 = [(L -1 k A)] \ ⌉ ⊙ F \ ⌉ , (34)
L k+1 = [(AU -1 k+1 )] ⌊ ⊙ F ⌊ . (35)
The filtering operation in (34) is equivalent to subtracting an upper triangular sparse matrix

S k from the matrix (L -1 k A) \ ⌉ , so it follows that U k+1 = [(L -1 k A)] \ ⌉ -S k .
Proceeding as in Proposition 1 we multiply this relation to the left by L k , so

L k U k+1 = L k [L -1 k A] \ ⌉ -L k S k
, and then apply the upper triangular filter to obtain, exploiting again Lemma 3,

[L k U k+1 ] \ ⌉ = [L k (L -1 k A) \ ⌉ ] \ ⌉ -[L k S k ] \ ⌉ = [A] \ ⌉ -[L k S k ] \ ⌉ . (36)
Similarly, equation (35

) yields L k+1 = [(AU -1 k+1 )] ⌊ ⊙F ⌊ = [(AU -1 k+1 )] ⌊ -T k ,
where T k represents the (sparse) lower triangular matrix of dropped entries. Then, proceeding in the same way as above, we get

L k+1 U k+1 = [(AU -1 k+1 )] ⌊ U k+1 -T k U k+1 and [L k+1 U k+1 ] ⌊ \ = [[(AU -1 k+1 )] ⌊ U k+1 ] ⌊ \ -[T k U k+1 ] ⌊ \ = A ⌊ \ -[T k U k+1 ] ⌊ \ . (37)
We have just shown the following two relations

[A -L k U k+1 ] \ ⌉ , = [L k S k ] \ ⌉ (38) [A -L k+1 U k+1 ] ⌊ \ . = [T k U k+1 ] ⌊ \ (39)
Though these relations do not establish convergence, they tell us something in the case when the process converges.

Proposition 2. Assume that a constant filter F is used and let S k and T k be the upper triangular and lower triangular matrices, respectively, of the dropped entries from U k+1 and L k+1 due to the filtering operation. Assume that the matrices L k and U k converge to L, Ū , respectively, and that Ū is nonsingular, Then,

(i) The sequence of matrices S k converges to S = [ L-1 A] \ ⌉ -Ū . (ii) The sequence of matrices T k converges to T = [A Ū -1 ] ⌊ -L.
(iii) The sequences of matrices

(A-L k U k+1 ) ⌊ \ and (A-L k U k ) \ ⌉ converge to [ LS] \ ⌉ and [T Ū ] ⌊ \ , respectively.
Proof. To prove (i) note that under the assumptions, Equation (38) implies that (L k S k ) \ ⌉ converges to [A -L Ū ] \ ⌉ . We now invoke Lemma 3 again and observe that:

[L -1 k (L k S k ) \ ⌉ ] \ ⌉ = [S k ] \ ⌉ = S k
The left-hand side, and therefore also S k , converges to the following limit

[ L-1 [A -L Ū ] \ ⌉ ] \ ⌉ = [ L-1 [A -L Ū ]] \ ⌉ = [ L-1 A] \ ⌉ -Ū .
The proof of (ii) is similar. The proof of (iii) is a direct consequence of relations (38-39).

The convergence of the matrices of dropped entries may seem counter intuitive. However, note that these matrices represent the cumulative matrix of dropped entries. Our experience indicates that the algorithm often converges rapidly to a certain limit. With a trivial filter equal to a diagonal matrix, the process generally converges in one or two steps. We ran a few tests with finite difference discretizations of a Laplacien on an n × n × n 3-D mesh, by taking the filter to be the original pattern of the matrix. We observed that the process converges in about N steps. In fact the norm L k -L k+1 F + U k -U k+1 F drops abruptly at around n steps.

Practical Implementations

Two key tools must be exploited in any implementation of the methods discussed so far. One concerns basic sparse-sparse computations, i.e., computations involving products of sparse matrices with sparse vectors or matrices. The other key tool concerns effective ways of solving sparse triangular systems with sparse right-hand sides. Basic sparse-sparse computations have been covered in detail in the literature, see e.g., [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF]. Techniques for solving sparse triangular systems with sparse right-hand sides are key ingredients of sparse direct solvers and so effective methods for this task have already been exploited. These are based on 'topological sorting', see, for example, [START_REF] Gilbert | Sparse partial pivoting in time proportional to artithmetic operations[END_REF][START_REF] Tinney | Sparse vector methods[END_REF][START_REF] Chan | Partial matrix refactorization[END_REF]. Thus, an LU factorization can be implemented as a sequence of sparse triangular solves where the right-hand sides are sparse. Because the right-hand side is sparse, the solution is also sparse, and it is easy to find a sequence in which the unknowns must be determined by examining a tree structure issued from the sparsity pattern of both the right-hand side and the matrix. This well-known idea is summarized in the next subsection, and details can be found in, e.g., [START_REF] Gilbert | Sparse partial pivoting in time proportional to artithmetic operations[END_REF][START_REF] Tinney | Sparse vector methods[END_REF][START_REF] Davis | Direct methods for sparse linear systems[END_REF]. What is less undertood is how to adapt this idea to the context of approximate triangular solves. This is taken up in Section 4.2

Background: Solving sparse-sparse triangular systems

A sparse-sparse triangular linear system of equations is a linear system whose coefficient matrix as well as the right-hand side are sparse. It is common to model the tasks involved in a sparse Figure 1. A small 10 × 10 sparse lower triangular matrix and its corresponding DAG triangular system with the help of a a Directed Acyclic Graph (DAG) which represents the binary relation: unknown i must be determined before unknown j. An illustration is given in Figure 1.

It is useful to briefly recall a well-known technique used in sparse matrix methods for solving sparse-sparse triangular systems. A standard column version of a lower triangular system, such as the one shown in Figure 1 would first compute the first entry, x 1 of the solution. Then a multiple of the first column is subtracted from the right-hand side. The second unknown is determined and the process continues until all components have been calculated.

Assume now that the right-hand side is sparse. In this case, the solution vector will remain sparse during these column subtractions so at each step we just need to consider the next nonzero entry in the running right-hand-side b. The question is how to determine the next entry each time. Consider the example of the figure in the situation when b = e 1 , the first column of the identity matrix. After x 1 is determined, the first column multiplied by x 1 is subtrated from b = e 1 . This introduces 'fill-ins' in positions 2 and 4 of b. The next nonzero entry is now b 2 and we determine x 2 . This introduces fill-ins in locations 5 and 6. At this point the next nonzero entry is b 4 , so we determine x 4 , etc.

Using matlab notation, the general algorithm is as follows (note that the first step copies the right-hand side onto the solution vector):

0. x := b, j = 0. 1. While (j < N ) Do 2.
Determine j := index of next nonzero entry of x 3.

Compute x j := x j /L jj 4.

Compute x(j + 1 : N ) := x(j + 1 : N ) -x j * L(j + 1 : N, j)

End

The problem addressed by what is referred to a topological sorting is step 2 of the algorithm, namely finding the order in which to determine the unknowns. A topological sorting of the entries is obtained from the DAG. If b = e 1 then we perform a Depth-First Search traversal of the graph (visit subtrees, then visit current node) starting from the root at node 1. The topological sort of the nodes is then just a reverse post-order traversal of the tree, see [START_REF] Davis | Direct methods for sparse linear systems[END_REF].

For the above example, the post order traversal starting from node 1, is 5, 6, 2, 9, 7, 4, 1 so the topological order is 1, 4, 7, 9, 2, 6, 5. Note that not all the graph is traversed and that we get a valid order that is different from the natural one discussed above. For a general situation, where the right-hand side has more than one nonzero entry we will proceed in the same way but we now need to visit the tree from several roots in turn. For example, if the right-hand side for the above example is the vector e 1 + e 3 , then we will do a traversal from node 1 and node 3. It is important to realize that the order is not important, i.e., starting the traversal from node 3 first and then node 1 or the other way around would both yield a valid order.

Approximate solutions of sparse-sparse triangular systems

So far we discussed only the exact solution case but we need to consider the situation of approximate solves since we are interested in approximate LU factorizations, not exact ones Numer. Linear Algebra Appl. 2009; 0:1-1 Prepared using nlaauth.cls as in sparse Gaussian Elimination. This particular situation does not seem to have been given much attention in the literature.

Suppose at first that b = e 1 again. In this situation we can just truncate the DFS to a certain depth level. This means that each node will get a level number inherited from its parent in the traversal, and when the level is higher than a certain maximum, denoted by 'MaxLev' thereafter, the traversal is curtailed at that point. In the above example, a traversal with a MaxLev of 2 would yield 5, 6, 2, 7, 4, 1. This is easy to justify with the common model of level-of fill ILU factorizations. For model problems, the higher the level the smaller the entry being considered, and so for example the element x 9 is likely to be smaller than its parent x 7 which in turn is smaller than its parent x 4 , etc.

A complication arises when we consider sparse right-hand sides with more than one nonzero entry. Consider again the situation when b = e 1 + e 3 in the above example. Then with respect to the traversal from node 1, x 7 has a level of 2, and x 9 a level of 3. With respect to the traversal from node 3, x 7 has a level of 1, and x 9 a level of 2. Which ones should we use? If we start the traversal from node 1 first, and if MaxLev=2, then 7 will be kept (lev=2) and 9 will be dropped (lev=3). Then when we will start the visit from node 3, the DFS traversal will see that node 7 was already visited and will ignore it.

We could resolve the issue by solving several systems with different right-hand sides, namely one for each nonzero entry, but this increases cost. A simpler solution, to prevent this is to order the traversal from higher to lower numbered nodes (in the example: first 3, then 1). This works simply because if j > i then i cannot be a descendent of j.

A simplified iteration

Consider the situation when A is perturbed as à = A + E where E is a very sparse matrix of small norm. Thus, we already have an approximate factorization LU = A + R, where R is the residual matrix in the approximation. After the perturbation we have à = LU + E + R. A natural idea suggested in [START_REF] Birken | Preconditioner updates applied to CFD model problems[END_REF] is to try to find a correction to L, U so that the residual norm remains the same as that of the original one. In other words we seek perturbations X L , X U to L, U respectively so that

(L + X L )(U + X U ) = (A + E) + R which leads to X L U + LX U + X L X U -E = 0. ( 40
)
Note that this is the same equation as (4) with R simply replaced by E. We could approximately solve these equations by performing one iteration of the corresponding variant of [START_REF] Golub | Matrix Computations[END_REF][START_REF] Grote | Parallel preconditioning and approximate inverses on the connection machine[END_REF]:

U k+1 = U k + [L -1 k E] \ ⌉ (41)
L k+1 = L k + [EU -1 k+1 ] ⌊ \ . (42)
Observe that when E is lower triangular, for example E = e i e T j with i > j, then the first correction, (41), is vacuous, i.e., U k+1 = U k , while the second correction gives

L k+1 = L k + [EU -1 k ] ⌊ \ . Approximating [EU -1 k ] ⌊ \ by [E(U ) -1
] ⌊ \ yields the approximation used in [START_REF] Birken | Preconditioner updates applied to CFD model problems[END_REF]. The tests reported in [START_REF] Birken | Preconditioner updates applied to CFD model problems[END_REF] with this approach were for a specific application in CFD. Our experiments show that for a general matrix and a general perturbation E then this very simple scheme does not improve the factorization to a satisfactory level.

An approach which we found to be more generally successful is one that is based on a simplification iterations [START_REF] Golub | Matrix Computations[END_REF][START_REF] Grote | Parallel preconditioning and approximate inverses on the connection machine[END_REF]. This simplification amounts to taking the maximum depth level MaxLev in the triangular solves to be zero, resulting essentially in replacing the triangular matrices by their diagonals. A few other minor changes are made to reduce cost. The algorithms is as follows where we assume that A = L 0 U 0 + R 0 is the initial factorization.

Algorithm 3 (Simplified ITALU)

1.

For k := 0, 1, . . . , Do: 2.

Compute

R k = A -L k U k 3.
Apply dropping to R k . 4.

Compute

U k+1 := U k + (R k ) \ ⌉ 5.
Compute

L k+1 := L k + (R k ) \ ⌉ * ((U k )) -1 6.
End

Note that Line 4 substitutes for

U k+1 = U k +(L -1 k R k ) \ ⌉ , i.e.
, it is the result of approximating L k by the identity matrix. Similarly, Line 5 results from approximating the current U k by its main diagonal. Note also that unlike Algorithm 2, the residual matrix R k is only computed once per iteration. This contrasts with our original algorithm where R k is recomputed after each update (L k or U k ). Also dropping is applied to R k and not to the corrections X U and X L .

If one iteration is performed, the cost of the above algorithm is quite modest. The computation of R k and the 'pruning' step (steps 1 and 2) are basically negligible because they only require keeping track of entries that are dropped during the original LU factorization. The only cost is that of the logic to test whether an entry is to be dropped and the additional overhead in building the matrix R k (memory allocation, pointer set-ups, etc). Steps 3 and 4 are equally inexpensive. They correspond to adding a hopefully very sparse matrix to the current factors L k and U k .

The remarkable aspect of the above Algorithm is that it works quite well in spite of its low cost, in the situation when the perturbation matrix is relatively small. A small number of iterations, e.g, 1 to 3, is usually sufficient to yield good improvements. An illustration is shown in the numerical experiments section, see Section 5.3.

Applications to linear and stationary problems

To illustrate the behavior of the incremental ILU, we consider in this section a few problems involving 3D-convection-diffusion like finite difference matrix. All experiments were performed in Matlab.

Alternating Lower-Upper corrections as a means to obtain an ILU

We will illustrate a method for extracting an ILU factorization by starting from the inaccurate SSOR defined below (or an ILU(0)) and then improving the factorization progressively by using one step of the procedure ITALU described in Section 3.

This experiment works on a stationary problem, specifically a matrix is generated from a convection-diffusion problem on a 15 × 15 × 10 regular grid. This yields a problem of size n = 2250. We use finite difference discretization and generate a problem with constant convection terms in each direction. So, in each direction, the discretized operator yields a tridiagonal matrix tridiag[-1 + α, 2, -1 -α].

In the first test α is taken equal to 0.1 and the matrix was shifted by 0.3 to make it indefinite. This shift introduces two negative eigenvalues and also two eigenvalues very close to zero as can seen from a run with MATLAB's eigs function:

λ 1 = -0.112843039478229, λ 2 = -0.000397969025052, λ 3 = 0.000397969025052, λ 4 = 0.113638977528334, λ 5 = 0.122450476821596 • • •
In addition the condition number estimated by Matlab was 1.68e+05. The right-hand side of the system to solve is generated artificially as b = A * ones(n, 1) in Matlab notation.

Figure 2 shows the performance of GMRES with 4 different preconditioners. The first preconditioner is SSOR with ω = 1. In fact this is simply the factorization A ≈ L 0 U 0 with L 0 = I + A ⌊ \ (A) -1 and U 0 = A \ ⌉ . The second preconditioner takes these two factors and generates an improved factorization A ≈ L 1 U 1 by performing one iteration of the ITALU algorithm. Dropping is performed in this procedure as follows. First a relative tolerance is set to droptol = 0.2, the effect of which is to ignore entries in any generated column that is less than droptol times the norm of the original row. In addition, as a precaution against excessive fill-in, only the largest lfil=10 entries are kept (however in this example, this second criterion does not seem to have been triggered). The results, in the form of residual norm versus iteration count, are shown on the left panel of Figure 2.

In the last run, corresponding to ITALU(3), the resulting preconditioner yielded a convergence behavior that is identical with that obtained from the previous one (i.e., from ITALU(2)). Therefore, we changed the parameters to show one more run, namely one corresponding to droptol = 0.1 and lfil = 20 (starting again from the preconditioner obtained in the previous run, i.e., the one labelled ITALU(2) in the figure). The numbers shown at the end of each curve (3.64, 1.97, 1.78, 1) are the "fill-factors" for each preconditioner, i.e., the ratios of the number of nonzero entries in L + U over that of A. These measure the amount of memory required to store the preconditioner and give a rough measure of the computational effort to construct it.

In a second test, we repeated exactly the same experiment with a problem of the same size but with different coefficients: we took α = 0.2 and then shifted the matrix by the diagonal -0.5 I to make the problem more indefinite. In this case SSOR has trouble converging. The results are shown on the right panel of Recall that the test corresponding to ITALU(3) used droptol=0.1 and lfil=20 instead of droptol=0.2 and lfil=10. The fill-factor is now substantially larger for this last test, but this results in a much faster convergence.

Descent-type methods for computing an ILU

The Frobenius-based method MERLU described in Section 2.2 will now be illustrated on the second example seen in the previous section, i.e., the Convection-Diffusion problem, with the parameter α = 0.2 for the convection term and the shift -0.5.

The following test is similar to the one described above. First GMRES( 30) is attempted with the SSOR preconditioner. This yields the first convergence curve of the left side of Figure 3. Then two improved preconditioners are generated by Algorithm 1. The same dropping parameters are applied in Lines 5 and 10 as in the previous example. We took droptol = 0.2 and lfil=10. As before the number shown at the end of each curve corresponds to the fill-factor for each case. The last curve shows the performance of an algorithm which is a modification of Algorithm 1 in which the dropping is done differently. In this variation (labeled MERLU1), the dropping in Lines 5 and 10 is removed. Instead, we drop directly from the factors U and L after each update. From a practical point of view this is more interesting because fill-in is better controlled. From a theoretical point of view, we loose the guarantee that the factors generated will have a diminishing residual norm A -LU F . The last curve in the figure shows this result with the preconditioner obtained from 4 steps of this procedure starting with the SSOR factors as the initial guess. The dropping parameters are again droptol = 0.2, lfil = 10.

We repeated the experiment with less restrictive fill-in parameter by changing the value of lfil to 5 instead of 10. This means that fewer entries will be kept each time. The results are shown on the right side of Figure 3. The alternative procedure MERLU1 converges slowly but requires very little fill-in (a fill ratio of 1.31). In fact the last few iterations of the iterative LU procedure to compute the LU factors are not very effective here because we keep too few entries at each iteration. As a result an almost identical convergence curve is obtained if the number of iterations is reduced from 4 to 2. To remedy this, we tested yet another variant which allowed more fill in after each outer iteration. Specifically, we simply increase lfil by 2 after each iterative LU computation. The resulting method is labeled MERLU2 in the plot. As can be seen this variant achieves good convergence relative to the fill-in allowed.

Another set of experiments performed to show an interesting application of the incremental approach is one in which the ILU factorization starts from a basic one such as an SSOR 30) for a convection-diffusion problem with the SSOR and a few levels of improvements using the Frobenius-based approach. The test matrix is the same as the second one of the previous section. The right plot shows tests with different dropping strategies.

preconditioner, and then is gradually improved at the occasion of each GMRES restart. We do not report these experiments here.

Use of Topological Sort and Simplified ITALU

His experiment focuses on the issue of efficient implementations, see Section 4.1. Though our current implementations are in matlab, it is possible to get an idea of the cost of the techniques by means of a rough operation count. The experimental setting is as follows. We generate a problem similar to the ones of the previous section taking this time n x = n y = 70, n z = 1, α = 0.05 and the resulting matrix is shifted by σ = -0.01. The matrix size is now n = 4, 900. We use the same artificial right-hand side as before, namely, b = A * ones(n, 1).

Six experiments are run on this example. First we solve the system with GMRES(30)-ILUT using matlab's luinc with a droptol of 0.1. The convergence is slow and the wanted reduction of the residual norm by a factor of 1.e-08 was not reached in 200 steps.

We then perturb the matrix as follows. We generate the following perturbation using matlab commands:

z = zeros(n,1); v = z; k = fix(n/20); k1 = fix(k/2); z(k:k:n) = (0.15/n)*(k:k:n); v(k1:k:n) = (0.15/n)*(k1:k:n); E = sparse(v) * sparse(z)';
This matrix E is added to A to produce the perturbed matrix B. This is a rank-one perturbation to A, and so there exist more effective ways to solve the new linear system, via the Sherman-Morisson formula [START_REF] Golub | Matrix Computations[END_REF] but we selected this form of perturbing A only for the sake of reproducibility. Perturbations with other randomly generated matrices yielded similar results. We now want to solve the linear systems Bx = b, with the same right-hand side b as before. The next test is to just try using the same preconditioner to solve the new linear system. This usually works reasonably well when the perturbation is small and / or of small rank. In this case however, the convergence is quite slow, substantially slower than with the original system. We then tested the performance of ITALU with the topological sorting algorithm discussed in Section 4.1. We used a Maximum Level of 3. The left side of Figure 4 shows the performance of GMRES [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] with the 4 tests based on topological sorting. These runs correspond to: 1) Original problem solved with GMRES(30) using the preconditioner ILUT(0.1) on the original problem. 2) Performance of the same preconditioner-accelerator pair on the perturbed problem.

3) Performance of the corrected preconditioner using MaxLev=3 in the topological sort. One iteration is made. 4) Performance of the corrected preconditioner using MaxLev=3 in the topological sort. Two correction iterations made.

The first two curves seen on the right side of the same figure show the same result as those on the right side and they are reproduced here to facilitate comparisons. Note that the small differences are caused by the fact that the initial guesses are random vectors. The curves 3) and 4) are replaced by their analogues using the simplified version of ITALU as represented by Algorithm 3. The numbers at the end of each curve indicate the fill-factors for each case, see Section 5.2.

While the matlab codes used in the first experiment is expensive relative to the C-based luinc it is clear that the simplified version should be very inexpensive. To give an idea, the time to construct the initial preconditioner for this example by luinc is 0.08 sec. The time to compute the first correction is 0.0493 and that for the second correction 0.08 (all seconds). It is to be stressed that the luinc code of matlab is optimized whereas the computations of the corrections only used simple-minded matlab functions (e.g., triu,..).

An application in an evolutive computational fluid dynamics problem

This section discusses an application of the methods proposed in this paper to the numerical simulation of the variable density incompressible Navier-Stokes system given on a domain Ω ⊂ R 2 . These equations model the motion of two immiscible fluids with varying density, evolving in a domain Ω. This is the case of the Rayleigh-Taylor instabilities or the water bubble falling in air -see [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]. These phenomena vary locally both in the physical domain and in time and the corresponding operators are time dependent.

Governing equations and discrete problems

The variable density incompressible Navier-Stokes system in Ω ⊂ R 2 can be written as:

∂ t ρ + div x (ρu) = 0, ( 43
)
ρ ∂ t u + (u • ∇ x )u + ∇ x p -µ∆ x u = f , (44)
div x u = 0. ( 45
)
Here ρ(t, x) ≥ 0 stands for the density of a viscous fluid whose velocity field is u(t, x) ∈ R 2 .

Copyright c 2009 John Wiley & Sons, Ltd.

Numer. Linear Algebra Appl. 2009; 0:1-1 Prepared using nlaauth.cls

The description of the external force is embodied into the right hand side f (t, x) of (44) and µ > 0 is the viscosity of the fluid. The unknowns depend both on time t ≥ 0 and position x ∈ Ω ⊂ R 2 . The third unknown of the problem is the pressure p(t, x) ∈ R; it can be seen as a Lagrange multiplier associated to the incompressibility constraint (45). The problem is well posed when the equations are completed by initial and boundary conditions.

In [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], different numerical methods have been proposed for the transport equation ( 43) and for evaluating the evolution of the velocity driven by (44-45). To be more specific, thanks to a time-splitting (the "Strang splitting" [START_REF] Strang | On the construction and comparison of difference schemes[END_REF]), it is possible to solve (43) for a given velocity by using a Finite Volume approach (see e.g. [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]) which is efficient when dealing with a pure convection equation and then, compute the divergence free solution of the momentum equation ( 44) by exploiting the advantages of Finite Element (FE) methods (see e.g. [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Ern | Éléments finis: théorie, applications, mise en oeuvre[END_REF]).

For this, the domain Ω is approximated by a computational domain Ω h , discretized by a conforming and isotropic set of triangles T h , with mesh-size h. The FE spaces V h ⊂ (H 1 0 (Ω h )) 2 for the velocity u h and Q h ⊂ L 2 0 (Ω h ) for the pressure p h must verify the "inf-sup" or LBB condition. In what follows we shall use the classical P 2 /P 1 elements in order to obtain a "saddle point problem" corresponding to the momentum equation (44-45). We also need to define a suitable discrete approximation of the (given) function ρ, compatible with the discretization of the velocity and the pressure. We use P 1 elements for the natural choice

ρ h ∈ Qh ⊂ L 2 (Ω h ).
Let us denote by ∆t the time step and t n = n ∆t, n ≥ 0. For the time discretization of (44-45), we use a classical backward difference method for the time derivative and a semi-implicit linearization scheme to treat the nonlinear convection term in the momentum equation (44). This scheme, already used in [START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF], has a second-order accuracy in time. Specifically, given the approximations u n , u n-1 and ρ * (ρ * = ρ n or ρ n+1 depending on the time step of the Strang splitting), at time t n+1 we compute (u n+1 , p n+1 ) by solving

ρ * 3u n+1 -4u n + u n-1 2∆t + (ū n+1 • ∇ x )u n+1 -µ∆ x u n+1 + ∇ x p n+1 = f n+1 , (46) div x u n+1 = 0. ( 47
)
Here, ūn+1 = 2u nu n-1 is the linear second-order extrapolation of the velocity field at t n+1 . The algebraic version of discrete variational form of (46-47) may be written in the following saddle point problem:

A n+1 B T B 0 U P = F 0 , (48)
where A n+1 is a nonsymmetric matrix associated with the convection-diffusion term, B is the matrix associated with the divergence operator, and the vectors U, P, F correspond, respectively, to the unknowns u n+1 h , p n+1 h and to the given function f n+1 h plus other explicit terms depending on the velocity. Here, A n+1 changes from one time step to the next, due to the evolution of the density and to the linearization of the convective term. Specifically, we have

A n+1 = 3 2∆t M + 1 Re L + D,
where M is the velocity mass matrix depending on ρ * h , L is the Laplacian matrix for the velocity and the matrix D corresponds to the convective term and depends on ρ * h and ūn+1 h . Typically, the matrices M and D change at each time step and are not symmetric. Finally, recall that the Reynolds number Re is proportional to the size of the domain, a reference density and a reference velocity and inversely proportional to the dynamic viscosity µ.

In the case of homogeneous density ρ(t, x) = ρ, ∀t > 0, the saddle point problem (48) has been studied in many papers, see e.g. [START_REF] Calgaro | A preconditioner for generalized saddle-point problems: Application to 3d stationary Navier-Stokes equations[END_REF] and references therein. In these references block diagonal or block triangular preconditioning strategies are advocated for the system matrix of (48). The eigenvalue analysis in these references shows that the eigenvalues of the preconditioned system are bounded independently of the mesh size h of the underlying grid T h and numerical tests show that the convergence rates for GMRES iterations deteriorates roughly like the Reynolds number Re.

The problem considered is to solve the consecutive systems (48) by taking advantage of earlier systems. In this section, we study the influence of the alternating lower-upper correction, Algorithm 2, as preconditioner of the nonsymmetric matrix A n+1 , in the case of a block triangular preconditioner

A B T 0 S or A 0 B S ,
applied to the saddle point problem, where A is an approximation of the nonsymmetric matrix A n+1 and S is an approximation of the pressure Schur complement BA -1 n+1 B T . Since FGMRES, the flexible variant of GMRES, allows variations in the preconditioner, it is used as the accelerator in our implementation. In FGMRES, preconditioning is applied to the right:

A n+1 B T B 0 A B T 0 S -1 W U W P = F 0 , with U P = A B T 0 S -1 W U W P .
We focus our study only on the preconditioner A. For the approximation of the pressure Schur complement BA -1 n+1 B T , we choose two classical preconditioners: the pressure-mass preconditioner [29] S -1 = M -1 p and the Cahouet-Chabard preconditioner [6] S -1 = 2∆t 3 M -1 p + ReL -1 p , where M p is the mass matrix for the pressure and L p is the discretization of the Laplace operator for the pressure.

The ITALU algorithm for saddle point problems

The efficiency of a preconditioner for the saddle point problem (48) can be tested numerically by verifying that the convergence rates for FGMRES iterations depend only mildly on the Reynolds number and, most importantly, that the number of iterations does not grow as the mesh size is reduced.

We want to evaluate the ability of the Algorithm 2 to recover analytical solution, checking the rates of convergence and comparing the mean number of iterations required by standard ILUT preconditioning (using Matlab luinc function with droptol = 10 -5 ) with those required by the ITALU algorithm. In particular, we use the same droptol and we increase the lfil and itlu parameters when h decreases. Here, itlu is the number of iterations performed by the ITALU algorithm. We check for each test the rates of convergence of the velocity and pressure considering the example given in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]:

       ρ ex (t,
Computations are performed for 0 ≤ t ≤ 0.5 on the square [-1, 1] 2 , using a structured mesh T h . The convergence results are displayed with respect to h and for different time steps ∆t. Table I shows the mesh size, the corresponding dimensions of vectors U and P and the number of nonzero entries in the matrix A n+1 , which is block-diagonal. For the operation with the matrix S -1 , we always use the L and U factors obtained from Matlab's ILUT (the luinc procedure), with the droptol set to droptol= 10 -3 and no pivoting.

The results in Figure 5 correspond to the pressure-mass preconditioner S -1 = M -1 p for the saddle point problem, when Re = 1 and h = h 1 , h 2 , h 3 . In order to compare the results of a standard approach and the new procedure presented in section 3, two algorithms are used for operations with the matrix A: the first one is the standard ILU preconditioner of Matlab luinc, by setting droptol = 10 -5 and no pivoting. The second one is the ITALU procedure with the same droptol. In this figure and the next two ones, the mean value of the fill-ratio of the preconditioner A is displayed for each case. Each of these plots shows iteration counts (vertical axis) versus -log 10 (∆t) (horizontal axis). The left plots show performances of the Matlab luinc(drop) preconditioning and the right plots those of the ITALU procedure.

It seems obvious that the parameters lfil and itlu should not be fixed. They should vary according to the mesh size h. In particular, we set lfil = 20,25,30 (or 40) and itlu = 3,4,5 respectively, when the mesh size h decreases. Smaller values of lfil or itlu, result in a greater number of FGMRES iterations on average. The results are shown on the right side of Figure 5. Near h = h 3 we indicate the lfil used. As expected, we can observe that the lfil parameter is very important to obtain a rate of convergence independent of the mesh size.

Using the same values for h, we test the Cahouet-Chabard preconditioner S -1 = 2∆t 3 M -1 p + ReL -1 p , for the saddle point problem, when Re = 1 and Re = 1000. This choice is necessary because for higher Reynolds numbers the performances of the pressure-mass preconditioner for the Schur complement are not satisfactory.

The results in Figure 6 (resp. Figure 7) correspond to Re = 1 (resp. Re = 1000). The results in the left side of these figures correspond to the standard ILUT preconditioner and those in the right side to the ITALU preconditioner. If Re = 1, we take the same parameters droptol, lfil and itlu as in the previous test case. As in the preceding case, the importance of the lfil parameter can be observed for small values of h. With higher Reynolds numbers, we must increase the parameters lfil and itlu in the ITALU preconditioner. In particular, we set lfil = 30,40,50 (or 70) and itlu = 3,4,5 respectively, when the mesh size h decreases. In this case, choosing higher lfil value does not reduce the mean number of iterations of FGMRES; the only consequence is an increased fill-in of A.

Overall, we stress that the fill-in of A is reduced in a significant way (up to 40%) when considering the ITALU procedure presented in Section 3 with respect to the classical approach, in particular when Re is very large. This crucial point will become very important when considering large problems for which memory costs will cause severe limitations.

However, to capture correctly the dynamics with higher Reynolds numbers, only small time steps ∆t must be chosen and in this case the performances of FGMRES are equivalent.

We also underline that, when we proceed through the sequence of problems, the number of iterations of FGMRES does not change very much with respect to the mean iteration number plotted in Figures 5, 6 and7 . Indeed, for both ITALU and luinc algorithms, the FGMRES iterations fluctuate between approximatively +10 % of the mean iteration number, at the beginning of the time stepping iterations, and decrease until -10% at the end of simulations.

Conclusion

Though much is known about finding effective preconditioners to solve general sparse linear systems which arise in real-life applications, little has been done so far to address the issue of updating such preconditioners. This paper proposed a few possible practical approaches. One option is to utilize techniques based on sparse approximate inverses. These methods can yield very inexpensive ways to correct LU factors. They rely on products of sparse matrices with sparse vectors which have been exploited in the literature, see, e.g., [START_REF]Approximate inverse preconditioners via sparse-sparse iterations[END_REF].

This paper also considered an alternating procedure, named ITALU, which corrects the factors L and U alternatingly. Experiments show that the method can perform quite well. The main ingredient in the construction of ITALU preconditioners is the solution of a sparse triangular system with a sparse right-hand side matrix. Computations of this type constitute the innermost kernels of sparse direct solvers, see, e.g., [START_REF] Davis | Direct methods for sparse linear systems[END_REF][START_REF] Gilbert | Sparse partial pivoting in time proportional to artithmetic operations[END_REF][START_REF] Tinney | Sparse vector methods[END_REF]. This paper has shown how to exploit the idea of levels in the graph traversal to adapt these techniques for the case when only an approximate solution is desired. At one extreme one can obtain a very inexpensive and surprisingly effective update of the ILU factorization by solving these systems in a very approximate manner. The corresponding technique, called the simplified ITALU, yielded nice improvements on the examples we tested. One can make the corrections as inexpensive as desired, by adjusting the depth levels in the sparse triangular solve. At the other extreme, one can obtain full-fledged methods for building accurate factorizations starting from, say, a diagonal or an SSOR factorization.

The possibility of computing updates of ILU factorizations inexpensively opens up many interesting avenues. Many applications in CFD lead naturally to small variations in the matrix. A feature that is common to many physical applications is that the varying parts or the 'rough' parts of the domain are usually limited in size. In examples involving the interaction of two fluids or a solid body and a fluid, the change in the coefficient matrix is localized around the physical interfaces. In this situation, the largest entries of the residual matrix R k will be localized too and the costs of the updates will be quite small. It is also possible to devise strategies whereby an ILU factorization is built with a certain uniform tolerance for the error A -LU F and then the factors are improved locally only where the domain is likely to cause difficulties. In this way one can improve the quality of the preconditioner only in areas where it is more likely that a higher accuracy is required. These are but a few of the topics which

Copyright c 2009

 2009 John Wiley & Sons, Ltd.

Figure 2 .

 2 The smallest eigenvalues as computed by the function eigs are now λ 1 = -0.2240921324 λ 2 = λ 3 = -.1125800958 λ 4 = -0.0010680592 λ 5 = .0076089054 λ 6 = λ 7 = 0.0685056254 λ 8 = .1191209420 ... and Matlab's condition number estimator yields condest(A) ≈ 4.03E + 05.

Figure 2 .

 2 Figure 2. Performance of GMRES(30) for a convection-diffusion problem with the SSOR and 3 levels of improvements by the ITALU algorithm. Left and right plot are for two different problems.

Copyright c 2009

 2009 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2009; 0:1-1 Prepared using nlaauth.cls

Figure 3 .

 3 Figure3. Performance of GMRES[START_REF] Strang | On the construction and comparison of difference schemes[END_REF] for a convection-diffusion problem with the SSOR and a few levels of improvements using the Frobenius-based approach. The test matrix is the same as the second one of the previous section. The right plot shows tests with different dropping strategies.

Copyright c 2009

 2009 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2009; 0:1-1 Prepared using nlaauth.cls

Figure 4 .

 4 Figure 4. Performance of GMRES(30) for a convection-diffusion problem with one or two levels of improvements of the ITALU algorithm (left side) and the simplified version described by Alg. 3.

Figure 6 .

 6 Figure 6. Mean number of iterations required by FGMRES(30) for the saddle point problem using the Cahouet-Chabard preconditioner for the Schur complement, for different values of h and Re = 1.

Figure 7 .

 7 Figure 7. Mean number of iterations required by FGMRES(30) for the saddle point problem using the Cahouet-Chabard preconditioner for the Schur complement, for different values of h and Re = 1000.

  x, y) = 2 + x cos(sin t) + y sin(sin t), Prepared using nlaauth.cls Figure 5. Mean number of iterations required by FGMRES(30) for the saddle point problem using pressure-mass preconditioner for the Schur complement, for different values of h and Re = 1.

	u ex (t, x, y) =	-y cos t x cos t	,	(49)
	p ex (t, x, y) = sin x sin y sin t.
	Copyright c 2009 John Wiley & Sons, Ltd.			Numer. Linear Algebra Appl. 2009; 0:1-1

The fields ρ ex (t, x, y) and u ex (t, x, y) satisfy the mass conservation equation (

43

) identically and u ex (t, x, y) is solenoidal. The momentum equation (

44

) is satisfied with the body force:

f ex (t, x, y) =

(y sin t -x cos 2 t)ρ ex (t, x, y) + cos x sin y sin t -(x sin t + y cos 2 t)ρ ex (t, x, y) + sin x cos y sin t .

Table I .

 I Dimensions of the problem with respect to the mesh size.

	h	size(U ) size(P ) nnz(A n+1 )
	h 1 = 0.0312	1832	289	20578
	h 2 = 0.0156	7938	1089	88162
	h 3 = 0.0078	32258	4225	364642

† Work initiated while the third author was visiting the department of Mathematics of the University of Lille. Work of the third author was supported by the US Department Of Energy under contract DE-FG-08ER25841 and by the Minnesota Supercomputer Institute.

are left for future studies.