Approximate Inverse Techniques for Block-Partitioned Matrices.
This paper proposes some preconditioning options when the system matrix is in block-partitioned form.  This form may arise naturally, for example, from the incompressible Navier--Stokes equations, or may be imposed after a domain decomposition reordering.  Approximate inverse techniques are used to generate sparse approximate solutions whenever these are needed in forming the preconditioner. The storage requirements for these preconditioners may be much less than for incomplete LU factorization (ILU) preconditioners for tough, large-scale computational fluid dynamics (CFD) problems. The numerical experiments show that these preconditioners can help solve difficult linear systems whose coefficient matrices are highly indefinite.
Introduction
Consider the block partitioning of a matrix A, in the form
(1)
where the blocking naturally occurs due the ordering of the equations and the variables.
Matrices of this form arise in many applications, such as in the incompressible Navier-Stokes
equations, where the scalar momentum equations and the continuity condition
form separate blocks of equations. In the 2-D case, this is a system of the form
A =B @
pC A =B @
f u
f pC A (2)
where u and v represent the velocity components, and p represents the pressure. Here, the
B submatrix is a convection-diffusion operator, the F submatrices are pressure gradient
operators, and the E submatrices are velocity divergence operators.
Traditional techniques such as the Uzawa algorithm have been used for these problems,
often because the linear systems that must be solved are much smaller, or because there
are zeros or small values on the diagonal of the fully-coupled system. These so-called
segregated approaches, however, suffer from slow convergence rates when compared to
aggregated, or fully-coupled solution techniques.
Another source of partitioned matrices of the form (1) is the class of domain decomposition
methods. In these methods the interior nodes of a subdomain are ordered
consecutively, subdomain after subdomain, followed by the interface nodes ordered at
the end. This ordering of the unknowns gives rise to matrices which have the following
structure: 0
Typically, the linear systems associated with the B matrix produced by this reordering
are easy to solve, being the result of restricting the original PDE problem into a set of
independent and similar PDE problems on much smaller meshes. One of the motivations
for this approach is parallelism. This approach ultimately requires solution methods for
the Schur complement S. There is a danger, however, that for general matrices, B may
be singular after the reordering.
Much work has been done on exploiting some form of blocking in conjunction with
preconditioning. In one of the earlier papers on the subject, Concus, Golub, and Meurant
introduce the idea of block preconditioning, designed for block-tridiagonal matrices
whose diagonal blocks are tridiagonal. The inverses of tridiagonal matrices encountered
in the approximations are themselves approximated by tridiagonal matrices, exploiting an
exact formula for the inverse of a tridiagonal matrix. This was later extended to the more
general case where the diagonal blocks are arbitrary [4, 17]. In many of these cases, the
incomplete block factorizations are developed for matrices arising from the discretization
of PDE's [2, 3, 7, 17, 19] and utilize approximate inverses when diagonal blocks need to
be inverted. More recently, Elman and Silvester [13] proposed a few techniques for the
specific case of the Stokes and Navier-Stokes problems. A number of variations of Block-
preconditioners have also been developed [1, 9]. In these techniques the off-block
diagonal terms are either neglected or an attempt is made to approximate their effect.
This paper explores some preconditioning options when the matrix is expressed in
block-partitioned form, either naturally or after some domain decomposition type re-
ordering. The iterative method acts on the fully-coupled system, but the preconditioning
has some similarity to segregated methods. This approach only requires preconditioning
or approximate solves with submatrices, where the submatrices correspond to any combination
of operators, such as reaction, diffusion, and convection. It is particularly advantageous
to use the block-partitioned form if we know enough about the submatrices to
apply specialized preconditioners, for example operator-splitting and semi-discretization,
as well as lower-order discretizations.
Block-partitioned techniques also require the sparse approximate solution to sparse
linear systems. These solutions need to be sparse because they form the rows or columns
of the preconditioner, or are used in further computations. Dense solutions here will
cause the construction or the application of the preconditioner to be too expensive. This
problem is ideally suited for sparse approximate inverse techniques. The approximate
solution to the sparse system is found by
using an iterative method implemented with sparse matrix-sparse vector and sparse
vector-sparse vector operations. The intermediate and final solutions are forced to be
sparse by numerically dropping elements in x with small magnitudes. If the right-hand-
side b and the initial guess for x are sparse, this is a very economical method for computing
a sparse approximate solution. We have used this technique to construct preconditioners
based on approximating the inverse of A directly [6].
This paper is organized as follows. In Section 2 we describe the sparse approximate
inverse algorithm and some techniques for finding sparse approximate solutions with the
Schur complement. Section 3 describes how block-partitioned factorizations may be used
as preconditioners. The most effective of these are the approximate block LU factorization
and the approximate block Gauss-Seidel preconditioner. Section 4 reports the results of
several numerical experiments, including the performance of the new preconditioners on
problems arising from the incompressible Navier-Stokes equations.
Sparse approximate inverses and their use
It is common when developing preconditioners based on block techniques to face the need
to compute an approximation to the inverse of a sparse matrix or an approximation to
columns of the f in which both B and f are sparse. This is particularly the
case for block preconditioners for block-tridiagonal matrices [7, 19]. For these algorithms
to be practical, they must provide approximations that are sparse.
A number of techniques have recently been developed to construct a sparse approximate
inverse of a matrix, to be used as a preconditioner [5, 6, 8, 10, 15, 17, 18]. Many
of these techniques approximate each row or column independently, focusing on (in the
column-oriented case) the individual minimizations
where e j is the j-th column of the identity matrix. Such a preconditioner is distinctly
easier than most existing preconditioners to construct and apply on a massively parallel
computer. Because they do not rely on matrix factorizations, these preconditioners often
are complementary to ILU preconditioners [6, 22].
Previous approaches select a sparsity pattern for x and then minimize (4) in a least
squares sense. In our approach, we minimize (4) with a method that reduces the residual
norm at each step, such as Minimal Residual or FGMRES [20], beginning with a sparse
initial guess. Sparsity is preserved by dropping elements in the search direction or current
solution at each step based on their magnitude or criteria related to the residual norm
reduction. The final number of nonzeros in each column is guaranteed to be not more
than the parameter lfil. In the case of FGMRES, the Krylov basis is also kept sparse
by dropping small elements. To keep the iterations economical, all computations are
performed with sparse matrix-sparse vector or sparse vector-sparse vector operations.
For our application here, we point out that the approximate inverse technique for
each column may be generalized to find a sparse approximate solution to the sparse linear
problem by minimizing
possibly with an existing preconditioner M for A.
2.1 Approximate inverse algorithm
We describe a modification of the technique reported in [6] that guarantees the reduction
of the residual norm at each minimal residual step. Starting with a sparse initial guess,
the fill-in is increased by one at each iteration. At the end of each iteration, it is possible
to use a second stage that exchanges entries in the solution with new entries if this causes
a reduction in the residual norm. Without the second stage, entries in the solution cannot
be annihilated once they have been introduced. For the problems in this paper, however,
this second stage has not been necessary.
In the first stage, the search direction d is derived by dropping entries from the residual
direction r. So that the sparsity pattern of the solution x is controlled, d is chosen to have
the same sparsity pattern as x, plus one new entry, the largest entry in absolute value.
Minimization is performed by choosing the steplength
(Ad; Ad)
and thus the residual norm for the new solution is guaranteed to be not more than the
previous residual norm. The solution and the residual is updated at the end of this
stage. If A is indefinite, the normal equations residual direction A T r may be used as the
search direction, or simply to determine the location of the new fill-in. It is interesting
to note that the largest entry in A T r gives the greatest residual norm reduction in a
one-dimensional minimization. This explains why a transpose initial guess for the approximate
inverse combined with self-preconditioning (preconditioning r with the current
approximate inverse) is so effective for some problems [6].
There are many possibilities for the second stage. We choose to drop one entry in
x and introduce one new entry in d if this causes a decrease in the residual norm. The
candidate for dropping is the smallest absolute nonzero entry in x. The candidate to be
added is the largest absolute entry in the previous search direction (at the beginning of
stage 1) not already included in d. The previous direction is used so that the candidate
may be determined in stage 1, and an additional search is not required. The steplength
fi is chosen by minimizing the new residual norm
where e i is the i-th coordinate vector, x s is the entry in x to be dropped at position s
while fi is the entry to be added at position l (largest), and we have generalized
the notation so that b is the right-hand-side vector, previously denoted m j . Let A j denote
the j-th column of A. Then the minimization gives
which just involves one sparse SAXPY since b \Gamma Ax is already available as r, and one sparse
dot-product, since we may scale the columns of A to have unit 2-norm. It is guaranteed
that s 6= l since l is chosen from among the entries not including s.
The preconditioned version of the algorithm for minimizing kb \Gamma Axk 2 with explicit
preconditioner M may be summarized as follows. A is assumed to be scaled so that its
columns all have unit 2-norm. The number of inner iterations is usually chosen to be lfil
or somewhat larger.
Algorithm 2.1 Approximate inverse algorithm
1. Starting with some initial guess x, r := b \Gamma Ax
2. For do
3. t := Mr
4. Choose d to be t with the same pattern as x;
one entry which is the
largest remaining entry in absolute value
5. q := Ad
6. ff := (r;q)
7. r := r \Gamma ffq
8. x
9. s := index of smallest nonzero in abs(x)
10. l := index of largest nonzero in abs(t \Gamma d)
11. fi := (r
12. ~ r := r
13. If k~rk ! krk then
14. Set x s := 0 and x l := fi
15. r := ~ r
16. End if
17. End do
2.2 Sparse solutions with the Schur complement
Sparse approximate solutions with the Schur complement are often
required in the preconditioning for block-partitioned matrices. We will briefly describe
three approaches in this section: (1) approximating S, (2) approximating S \Gamma1 , and (3)
exploiting a partial approximate inverse of A.
2.2.1 Approximating S
To approximate S with a sparse matrix, we can use
~
where Y is computed by the approximate inverse technique, possibly preconditioned with
whatever we are using to solve with B. Since Y is sparse, ~
S computed this way is also
sparse. Moreover, since S is usually relatively dense, solving with ~
S is an economical
approach. Typically, a zero initial guess is used for Y . We remark that it is usually too
expensive to form Y by solving B approximately and then dropping small elements,
since it is rather costly to search for elements to drop. We also note that we can generate ~
column-by-column, and if necessary, compute a factorization of ~
S on a column-by-column
basis as well. The linear systems with ~ S can be solved in any fashion, including with an
iterative process with or without preconditioning.
2.2.2 Approximating S \Gamma1
Another method is to compute an approximation to S \Gamma1 using the idea of induced pre-
conditioning. is the (2,2) block of
\GammaS
we can compute a sparse approximation to it by using the approximate inverse technique
applied to the last block-column of A and then throwing away the upper block. In practice,
the upper part of each column may be discarded before computing the next column. In
our experiments, since the approximate inverse algorithm is applied to A, an indefinite
matrix in most of the problems, the normal equations search direction A T r is used in the
algorithm, with a scaled identity initial guess for the inverse.
2.2.3 Partial approximate inverse
A drawback of the above approach is that the top submatrix of the last block-column
is discarded, and that the resulting approximation of S \Gamma1 may actually contain very few
nonzeros. A related technique is to compute the partial approximate inverse of A in the
last block-row. This technique does not give an approximation to S \Gamma1 , but defines a
simple preconditioning method itself. Writing the inverse of A in the form,
we can then get an approximate solution to A
y
with
f
It is not necessary to solve accurately with B. Again, the normal equations search direction
is used for the approximate inverse algorithm in the numerical experiments. Some
results of this relatively inexpensive method will be given in Section 4.
Block-partitioned factorizations of A
We consider a sparse linear system
which is put in the block form,
y
For now the only condition we require on this partitioning is that B be nonsingular. We
use extensively the following block LU factorization of A,
I
in which S is the Schur complement,
As is well-known, we can solve (12) by solving the reduced system,
to compute y, and then back-substitute in the first block-row of the system (11) to obtain
x, i.e., compute x by
The above block structure can be exploited in several different ways to define preconditioners
for A. Thus, the block preconditioners to be defined in this section combine
one of the preconditioners for S seen in Section 2.2 and a choice of a block factorization.
Next, we describe a few such options.
3.1 Solving the preconditioned reduced system
A method that is often used is to solve the reduced system (14), possibly with the help
of a certain preconditioner M S for the Schur complement matrix S. Although this does
not involve any of the block factorizations discussed above, it is indirectly related to it
and to other well-known algorithms. For example, the Uzawa method which is typically
formulated on the full system, can be viewed as a Richardson (or fixed point) iteration
applied to the reduced system. The matrix S need not be computed explicitly; instead,
one can perform the matrix-vector product with the matrix S, via the following
sequence of operations:
1. Compute
2. Solve
3. Compute
If we wish to use a Krylov subspace technique such as GMRES on the preconditioned
reduced system, we need to solve the systems in Step 2, exactly, i.e., by a direct solver or
an iterative solver requiring a high accuracy. This is because the S matrix is the coefficient
matrix of the system to be solved, and it must be constant throughout the GMRES
iteration. We have experimented with this approach and found that this is a serious
limitation. Convergence is reached in a number of steps which is typically comparable
with that obtained with methods based on the full matrix. However, each step costs much
more, unless a direct solution technique is used, in which case the initial LU factorization
may be very expensive. Alternatively, a highly accurate ILU factorization can be employed
for B, to reduce the cost of the many systems that must be solved with it in the successive
outer steps.
3.2 Approximate block diagonal preconditioner
One of the simplest block preconditioners for a matrix A partitioned as in (1) is the
block-diagonal matrix
in which MC is some preconditioning for the matrix C. If as is the case for the
incompressible Navier-Stokes equations, then we can define I for example. An
interesting particular case is when C is nonsingular and MC = C. This corresponds to a
block-Jacobi iteration. In this case, we have
the eigenvalues of which are the square roots of the eigenvalues of the matrix C
Convergence will be fast if all these eigenvalues are small.
3.3 Approximate block LU factorization
The block factorization (12) suggests using preconditioners based on the block LU factor-
ization
in which
and
to precondition A. Here M S is some preconditioner to the Schur complement matrix S. If
we had a sparse approximation ~
S to the Schur complement S we could compute a preconditioning
matrix M S to ~
S, for example, in the form of an approximate LU factorization.
We must point out here that any preconditioner for S will induce a preconditioner for A.
As was discussed in Section 3.1 a notable disadvantage of an approach based on solving
the reduced system (14) by an iterative process is that the action of S on a vector must
be computed very accurately in the Krylov acceleration part. In an approach based on
the larger system (11) this is not necessary. In fact any iterative process can be used for
solving with M S and B provided we use a flexible variant of GMRES such as FGMRES
[20].
Systems involving B may be solved in many ways, depending on their difficulty and
what we know about B. If B is known to be well-conditioned, then triangular solves with
incomplete LU factors may be sufficient. For more difficult B matrices, the incomplete
factors may be used as a preconditioner for an inner iterative process for B. Further, if
the incomplete factors are unstable (see Section 4.2), an approximate inverse for B may
be used, either directly or as a preconditioner. If B is an operator, an approximation to it
may be used; its factors may again be used either directly or as a preconditioner. This kind
of flexibility is typical of what is available for using iterative methods on block-partitioned
matrices.
An important observation is that if we solve exactly with B then the error in this
block ILU factorization lies entirely in the (2,2) block since,
One can raise the question as to whether this approach is any better than one based on
solving the reduced system (14) preconditioned with M S . It is known that in fact the
two approaches are mathematically equivalent if we start with the proper initial guesses.
Specifically, the initial guess should make the x-part of the residual vector equal to 0 for
the original system (11), i.e., the initial guess is
with
This result, due to Eisenstat and reported in [16], immediately follows from (16) which
shows that the preconditioned matrix has the particular form,
Thus, if the initial residual has its x-component equal to zero then all iterates will be
vectors with y components only, and a GMRES iteration on the system will reduce to a
GMRES iteration with the matrix M \Gamma1
involving only the y variable.
There are many possible options for choosing the matrix M S . Among these we consider
the following ones.
ffl no preconditioning on S.
ffl precondition with the C matrix if it is nonsingular. Alternatively we can
precondition with an ILU factorization of C.
construct a sparse approximation to S and use it as a preconditioner. In
general, we only need to approximate the action of S on a vector, for example, with
the methods described in Sections 2.2.1 and 2.2.2.
The following algorithm applies one preconditioning step to
to get
y
Algorithm 3.1 Approximate block LU preconditioning
1. x
2. y
3. x :=
We have experimented with a number of options for solving systems with M S in step 2
of the algorithm above. For example, M S may be approximated with ~
computed by the approximate inverse technique. If this approximation is
used, it is possible to also use Y in place of B \Gamma1 F in step 3.
3.4 Approximate block Gauss-Seidel
By ignoring the U factor of the approximate block LU factorization, we are led to a form
of block Gauss-Seidel preconditioning, defined by
The same remarks on the ways to solve systems with B and ways to define the preconditioning
matrix M S apply here. The algorithm for this preconditioner is the same as
Algorithm 3.1 without step 3.
To analyze the preconditioner, we start by observing that
showing that the only difference with the preconditioned matrix (17) is the additional
in the (1,2) position. The iterates associated with the block form and those
of the associated Schur complement approach M \Gamma1
are no longer simply related.
However there are a few connections between (17) and (19). First, the spectra of the two
matrices are identical. This does not mean, however, that the two matrices will require
the same number of iterations to converge in general.
Consider a GMRES iteration to solve the preconditioned system M
Here, we take an initial guess of the form
in which x 0 is arbitrary. With this we denote the preconditioned initial residual by
Then GMRES will find a vector u of the form belonging to the Krylov
subspace
which will minimize kM . For an arbitrary u in the affine space
the preconditioned residual is of the form
and by (19) this becomes,
As a result,
ks
Note that ks
represents the preconditioned residual norm for the reduced
system for the y obtained from the approximation of the large system. We have
ks
which implies that if the residual for the bigger system is less than ffl, then the residual
obtained by using a full GMRES on the associated preconditioned reduced system
will also be less than ffl. We observe in passing that the second term in
the right-hand-side of (21) can always be reduced to zero by a post-processing step which
consists of forcing the first part of the residual to be zero by changing ffi (only) into:
Equivalently, once the current pair x; y is obtained, x can be recomputed by satisfying
the first block equation, i.e.,
This post-processing step requires only one additional B solve.
Assume now that we know something about the residual vector associated with m
steps of GMRES applied to the preconditioned reduced system. Can we say something
about the residual norm associated with the preconditioned unreduced system? We begin
by establishing a simple lemma.
Lemma 3.1 Let
Then, the following equality holds
Proof. First, it is easy to prove that
in which Y We now multiply both members of the above equality
I \Gamma Z to obtain,
'We now state the main result concerning the comparison between the two approaches.
Theorem 3.1 Assume that the reduced system (14) is solved with GMRES using the
preconditioner M S starting with an arbitrary initial guess y 0 and let s
the preconditioned residual obtained at the m-th step. Then the preconditioned residual
vector r m+1 obtained at the (m 1)-st step of GMRES for solving the block system (11)
preconditioned with the matrix M of (18) and with an initial guess u
in which x 0
is arbitrary satisfies the inequality
In particular if s
Proof. The preconditioned matrix for the unreduced system is of the form (22) with
S. The residual vector s m of the m-th GMRES approximation
associated with the reduced system is of the form,
in which ae m is the m-th residual polynomial, which minimizes kp(G)s 0 k 2 among all polynomials
p of degree m satisfying the constraint:
Consider the polynomial of degree m+ 1 defined by
It is clear that
The residual of um+1 , the m+1-st approximate solution obtained by the GMRES algorithm
for solving the preconditioned unreduced system minimizes p(Z)r 0 over all polynomials p
of degree m+ 1 which are consistent, i.e., such that Therefore,
Using the equality established in the lemma, we now observe that
The first matrix in the right-hand-side of the last equality is nothing but I \Gamma Z. Hence,
the residual vector r m+1 is such that
which completes the proof. 2
It is also interesting to relate the convergence of this algorithm to that of the block-diagonal
approach in the particular case when M case corresponds to a block
Gauss-Seidel iteration. We can exploit Young and Frankel's theory for 2-cyclic matrices
to compare the convergence rates of this and the block Jacobi approach. Indeed, in this
case, we have from (19) that
Therefore, the eigenvalues of this matrix are the squares of those of matrix I
associated with the block-Jacobi preconditioner of Section 3.2.
4 Numerical Experiments
This section is organized as follows. In Section 4.1 we describe the test problems and
list the methods that we use. In Section 4.2, we illustrate for comparison purposes the
difficulty of incomplete LU factorizations for solving these problems in a fully-coupled
manner. In Section 4.3, we make some comments in regard to domain decomposition
types of reorderings. In Section 4.4 we show some results of the new preconditioners on a
simple PDE problem. Finally, in Sections 4.5 and 4.6, we present the results of the new
preconditioners on more realistic problems arising from the incompressible Navier-Stokes
equations.
Linear systems were constructed so that the solution is a vector of all ones. A zero
initial guess for right-preconditioned FGMRES [20] restarted every 20 iterations was used
to solve the systems. The Tables show the number of iterations required to reduce the
residual norm by 10 \Gamma7 . The iterations were stopped when 300 matrix-vector multiplications
were reached, indicated by a dagger (y). The codes were written in FORTRAN 77
using many routines from SPARSKIT [23], and run in single precision on a Cray C90
supercomputer.
4.1 Test problems and methods
The first set of test problems is a finite difference Laplace equation with Dirichlet boundary
conditions. Three different sized grids were used. The matrices were reordered using a
domain decomposition reordering with 4 subdomains. In the following tables, n is the
order of the matrix, nnz is the number of nonzero entries, nB is the order of the B
submatrix, and nC is the order of the C submatrix.
The second set of test matrices were extracted from the example incompressible Navier-Stokes
problems in the FIDAP [14] package. All problems with zero C submatrix were
tested. In the case of transient problems, the matrices are the Jacobians when the Newton
iterations had converged. The matrices are reordered so that the continuity equations are
Grid n nnz nB nC
by
48 by 48 2209 10857 2116 93
by 64 3969 19593 3844 125

Table

1: Laplacian test problems.
ordered last. The scaling of many of the matrices are poor, since each matrix contains
different types of equations. Thus, we scale each row to have unit 2-norm, and then scale
each column the same way. The problems are all originally nonsymmetric except 4, 12,
14 and 32.
Matrix n nnz nB nC
Hamel flow
EX12 3973 79078 2839 1134 Stokes flow
Surface disturbance attenuation
EX23 1409 42761 1008 401 Fountain flow
coating
EX26 2163 74465 1706 457 Driven thermal convection
EX28 2603 77031 1853 750 Two merging liquids
species deposition
Radiation heat transfer
EX36 3079 53099 2575 504 Chemical vapor deposition

Table

2: FIDAP example matrices.
The third set of test problems is from a finite-element discretization of the square lid-
driven cavity problem. Rectangular elements were used, with biquadratic basis functions
for velocities, and linear discontinuous basis functions for pressure. We will show our
results for problems with Reynolds number 0, 500, and 1000. All matrices arise from a
mesh of 20 by 20 elements, leading to matrices of size having nnz =138,187
nonzero entries. These matrices have 3363 velocity unknowns, and 1199 pressure un-
knowns. The matrices are scaled the same way as for the FIDAP matrices-the problems
are otherwise very difficult to solve.
We will use the following names to denote the methods that we tested.
ILUT(nfil) and ILUTP(nfil) Incomplete LU factorization with threshold of nfil nonzeros
per row in each of the L and U factors. This preconditioner will be described in
Section 4.2.
PAR(lfil) Partial approximate inverse preconditioner described in Section 2.2.3, using
lfil nonzeros per row in M 2 .
ABJ Approximate block-Jacobi preconditioner described in Section 3.2. This preconditioner
only applies when C 6= 0.
ABLU(lfil) Approximate block LU factorization preconditioner described in Section 3.3.
The approximation (6) to S with lfil nonzeros per column of Y was used.
ABLU y(lfil) Same as above, but using Y whenever B needs to be applied in step
3 of Algorithm 3.1.
ABLU s(lfil ) Approximate block LU factorization preconditioner, using (7) to approximate
S \Gamma1 with lfil nonzeros per column when approximating the last block column
of the inverse of A.
ABGS(lfil) Approximate block Gauss-Seidel preconditioner described in Section 3.4.
The approximation (6) to S with lfil nonzeros per column of Y was used.
The storage requirements for each preconditioner are given in Table 3. The ILUT
preconditioner to be described in the next subsection requires considerably more storage
than the approximate block-partitioned factorizations, since its storage depends on n
rather than nC . Because the approximation to S \Gamma1 discards the upper block, the storage
for it is less than lfil \Thetan C . The storage required for ~
S is more difficult to estimate since it
is at least the product of two sparse matrices. It is generally less than 2 \Theta lfil \Theta nC ; Table
11 in Section 4.5 gives the exact number of nonzeros in ~
S for the FIDAP problems.
4.2 ILU for the fully-coupled system
We wish to compare our new preconditioners with the most general, and in our experi-
ence, one of the most effective general-purpose preconditioners for solving the fully-coupled
system. In particular, we show results for ILUT, a dual-threshold, incomplete LU factorization
preconditioner based on a drop-tolerance and the maximum number of new fill-in
elements allowed per row in each L and U factor. This latter threshold allows the storage
for the preconditioner to be known beforehand. Drop-tolerance ILU rather than level-fill
ILU is often more effective for indefinite problems where numerical values play a much
more important role. A variant that performs column pivoting, called ILUTP, is even
more suitable for highly indefinite problems.
Matrices Matrix locations
ABJ none none
S less than 2 \Theta lfil \Theta nC
ABLU y(lfil) ~
lfil \Theta nC
S less than 2 \Theta lfil \Theta nC

Table

3: Storage requirements for each preconditioner.
We use a small modification that we have found to often perform better and rarely
worse on matrices that have a wide ranging number of elements per row or column. This
arises for various reasons, including the fact that the matrix contains the discretization of
different equations. Instead of counting the number of new fill-ins, we keep the nonzeros
in each row of L and U fixed at nfil , regardless of the number of original nonzeros in that
row. We also found better performance when keeping nfil constant rather than having it
increase or decrease as the factorization progresses.
If A is highly indefinite or has large nonsymmetric parts, an ILU factorization often
produces unstable L and U factors, i.e., k(LU) \Gamma1 k can be extremely large, caused by the
long recurrences in the forward and backward triangular solves [11]. To illustrate this
point, we computed for a number of factorizations the rough lower bound
where e is a vector of all ones. For the FIDAP example matrix EX07 modeling natural
convection with order 1633 and 46626 nonzeros, we see in Table 4 that the norm bound
increases dramatically as nfil is decreased in the incomplete factorization. GMRES could
not solve the linear systems with these factorizations as the preconditioner. This matrix
we chose is a striking example because it can be solved without preconditioning.
log

Table

4: Estimate of k(LU) \Gamma1 k1 from ILUT factors for EX07.
To illustrate the difficulty of solving the FIDAP problems with ILUTP, we progressively
allowed more fill-in until the problem could be solved, incrementing nfil in multiples
of 10, with no drop tolerance. The results are shown in Table 5. For these types of prob-
lems, it is typical that very large amounts of fill-in must be used for the factorizations to
be successful. An iterative solution was not attempted if the LU condition lower bound
was greater than 10 . If a zero pivot must be used, ILUT and ILUTP attempt to complete
the factorization by using a small value proportional to the norm of the row. The
matrices were taken in their original banded ordering, where the degrees of freedom of a
node or element are numbered together. As discussed in the next subsection, this type of
ordering having low bandwidth is often essential for an ILU-type preconditioning-many
problems including these cannot be solved otherwise.
Matrix nfil
EX06 50

Table

5: nfil required to solve FIDAP problems with ILUTP.
We should note that ILUTP is occasionally worse than ILUT. This can be alleviated
somewhat by using a low value of mbloc, a parameter in ILUTP that determines how
far to search for a pivot. In summary, indefinite problems such as these arising from the
incompressible Navier-Stokes equations may be very tough for ILU-type preconditioners.
4.3 Domain decomposition reordering considerations
Graph partitioners subdivide a domain into a number of pieces and can be used to give the
domain decomposition reordering described in Section 1. This is a technique to impose
a block-partitioned structure on the matrix, and adapts it for parallel processing, since
B is now a block-diagonal matrix. This technique is also useful if B is highly indefinite
and produces an unstable LU factorization; by limiting the size of the factorization, the
instability cannot grow beyond a point for which the factorization is not useful. For
general, nonsymmetric matrices, the partitioner may be applied to a symmetrized graph.
In

Table

6 we show some results of ILUT(40) on the Driven cavity problem with
different matrix reorderings. We used the original unblocked ordering where the degrees
of freedom of the elements are ordered together, the blocked ordering where the continuity
equations are ordered last, and a domain decomposition reordering found using a simple
automatic recursive dissection procedure with four subdomains. This latter ordering found
nodes internal to the subdomains, and 882 interface nodes.
Re. Unblocked Blocked DD ordered
1000 78 y 51

Table

Effect of ordering on ILUT for Cavity problems.
The poorer quality of the incomplete factorization for the Driven cavity problems in
block-partitioned form is due to the poor ordering rather than instability of the L and U
factors; in fact, zero pivots are not encountered. For the problem with Reynolds number
0, the unblocked format produces 745,187 nonzeros in the strictly lower-triangular part
during the incomplete factorization (which is then dropped down to less than n\Theta
nonzeros) while the block-partitioned format produces 2,195,688 nonzeros, almost
three times more.
The factorization for the domain decomposition reordered matrices encounters many
zero pivots when it reaches the (2,2) block. These latter orderings do not necessarily
cause ILUT to fill-in zeros on the diagonal. Nevertheless, the substitution of a small pivot
described above seems to be effective here. The domain decomposition reordering also
reduces the amount of fill-in because of the shape of the matrix (a downward pointing
arrow). Combined with its tendency to limit the growth of instability, the results show
this reordering is advantageous even on serial computers.
In

Table

7 we compare the difficulty of solving the B and ~
S subsystems for the blocked
and domain decomposition reorderings of the Driven cavity problems. ~
S was computed
as ~
computed using the approximate inverse technique with lfil
of 30. Here we used ILUT(30) and only solved the linear systems to a tolerance of 10 \Gamma5 .
Solves with these submatrices in the block-partitioned preconditioners usually need to be
much less accurate. In most of the experiments that follow, we used unpreconditioned
iterations to a tolerance of 10 \Gamma1 or 100 matrix-vector multiplications to solve with B and
~
S. Other methods would be necessary depending on the difficulty of the problems. The
table gives an idea of how difficult it is to solve with B and ~
S, and again shows the
advantage of using domain decomposition reorderings for hard problems.
Re. Blocked DD ordered
1000 y y 7

Table

7: Solving with B and ~
S for different orderings of A.
4.4 Test results for the Laplacian problem
In

Tables

8 and 9 we present the results for the Laplacian problem with three different grid
sizes, using no preconditioning, approximate block diagonal, partial approximate inverse,
approximate block LU, and approximate block Gauss-Seidel preconditioners. Note that
in

Table

9, an lfil of zero for the approximate block LU and Gauss-Seidel preconditioners
respectively indicate the preconditioners
I
and
Grid NOPRE ABJ PAR
by
48 by 48 367 50 29 21 19 17
by 64 532 57 36 33 25 20

Table

8: Test results for the Laplacian problem.
Grid ABLU ABGS
by
48 by 48 17
by 64 19 20

Table

9: Test results for the Laplacian problem.
4.5 Test results for the FIDAP problems
For the block-partitioned factorization preconditioners, unpreconditioned GMRES, restarted
every 20 iterations, was used to approximately solve the inner systems involving B and ~
by reducing the initial residual norm by a factor of 0.1, or using up to 100 matrix-vector
multiplications. Solves with the matrix B are usually not too difficult because for most
problems, it is positive definite. A zero initial guess for these solves was used. The results
for a number of the preconditioners with various options are shown in Table 10. The best
preconditioner appears to be ABLU using Y for B \Gamma1 F is better than solving a system
with B very inaccurately. The number of nonzeros in ~
S is small, as illustrated by Table
11 for two values of lfil.
Matrix

Table

10: Test results for the FIDAP problems.
4.6 Test results for the Driven cavity problems
The driven cavity problems are much more challenging because the B block is no longer
positive definite, and in fact, acquires larger and larger negative eigenvalues as the Reynolds
number increases. For these problems, the unpreconditioned GMRES iterations with B
were done to a tolerance of 10 \Gamma3 or a maximum of 100 matrix-vector multiplications.
Again, ABLU y appears to be the best preconditioner. The results are shown in Table
lfil
EX26 13395 21468
EX36 13621 21063

Table

11: Number of nonzeros in ~
S.
ABLU ABLU y ABGS
1000 y y 164 118 y y

Table

12: Test results for the Driven cavity problems.
Conclusions
We have presented a few preconditioners which are defined by combining two ingredients:
(1) a sparse approximate inverse technique for obtaining a preconditioner for the Schur
complement or a part of the inverse of A, and (2) a block factorization for the full sys-
tem. The Schur complement S which appears in the block factorization is approximated
by its preconditioner. Approximate inverse techniques [6] are used in different ways to
approximate either S directly or a part of A \Gamma1 .
As can be seen by comparing Tables 5 and 10, we can solve more problems with the
block approach than with a standard ILU factorization. In addition, this is typically
achieved with a far smaller memory requirement than ILUT or a direct solver. The better
robustness of these methods is due to the fact that solves are only performed for small
matrices. In effect, we are implicitly using the power of the divide-and-conquer strategy
which is characteristic of domain decomposition methods. The smaller matrices obtained
from the block partitioning can be preconditioned with a standard ILUT approach. The
larger matrices use a block-ILU, and the glue between the two is the preconditioning of
the Schur complement.

Acknowledgements

The authors wish to acknowledge the support of the Minnesota
Supercomputer Institute which provided the computer facilities and an excellent environment
to conduct this research.


--R


Incomplete block matrix factorization preconditioning methods.
Iterative Solution Methods.
On approximate factorization methods for block matrices suitable for vector and parallel processors.
Iterative solution of large sparse linear systems arising in certain multidimensional approximation problems.
Approximate inverse preconditioners for general sparse matri- ces
Block preconditioning for the conjugate gradient method.
Approximate inverse preconditioning for sparse linear systems.
Parallelizable block diagonal preconditioners for the compressible Navier-Stokes equations
A new approach to parallel preconditioning with sparse approximate inverses.
A stability analysis of incomplete LU factorizations.
Multigrid and Krylov subspace methods for the discrete Stokes equa- tions
Fast nonsymmetric iterations and preconditioning for Navier-Stokes equations
FIDAP: Examples Manual
Parallel preconditioning and approximate inverses on the Connection Machine.
A comparison of domain decomposition techniques for elliptic partial differential equations and their parallel implementation.
On a family of two-level preconditionings of the incomplete block factorization type
Factorized sparse approximate inverse precon- ditionings I
Incomplete block factorizations as preconditioners for sparse SPD matrices.
A flexible inner-outer preconditioned GMRES algorithm
ILUT: A dual threshold incomplete LU factorization.
Preconditioned Krylov subspace methods for CFD applications
SPARSKIT: a basic tool kit for sparse matrix computations
--TR

--CTR
N. Guessous , O. Souhar, Multilevel block ILU preconditioner for sparse nonsymmetric M-matrices, Journal of Computational and Applied Mathematics, v.162 n.1, p.231-246, 1 January 2004
Kai Wang , Jun Zhang, Multigrid treatment and robustness enhancement for factored sparse approximate inverse preconditioning, Applied Numerical Mathematics, v.43 n.4, p.483-500, December 2002
Prasanth B. Nair , Arindam Choudhury , Andy J. Keane, Some greedy learning algorithms for sparse regression and classification with mercer kernels, The Journal of Machine Learning Research, 3, 3/1/2003
Edmond Chow , Michael A. Heroux, An object-oriented framework for block preconditioning, ACM Transactions on Mathematical Software (TOMS), v.24 n.2, p.159-183, June 1998
Howard C. Elman , Victoria E. Howle , John N. Shadid , Ray S. Tuminaro, A parallel block multi-level preconditioner for the 3D incompressible Navier--Stokes equations, Journal of Computational Physics, v.187 n.2, p.504-523, 20 May
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
sparse approximate inverse;block-partitioned matrix;navier-stokes;schur complement;preconditioning