Dynamically Relaxed Block Incomplete Factorizations for Solving Two- and Three-Dimensional Problems.
To efficiently solve second-order discrete elliptic PDEs by Krylov subspace-like methods, one needs to use some robust preconditioning techniques. Relaxed incomplete factorizations (RILU) are powerful candidates.\ Unfortunately, their efficiency critically depends on the choice of the relaxation parameter $\omega$ whose "optimal" value is not only hard to estimate but also strongly varies from one problem to another. These methods interpolate between the popular ILU and its modified variant (MILU). Concerning the pointwise schemes, a new variant of RILU that dynamically computes variable $\omega=\omega_i$ has been introduced recently. Like its ancestor RILU and unlike standard methods, it is robust with respect to both existence and performance. On top of that, it breaks the problem-dependence of ``$\omega_{opt}$."\ A block version of this dynamically relaxed method is proposed and compared with classical pointwise and blockwise methods as well as with some existing "dynamic" variants, showing that with the new blockwise preconditioning technique, anisotropies are handled more effectively.
Introduction
. As model problem, we take the following self-adjoint second
order elliptic PDE
in\Omega with suitable boundary conditions on
where\Omega denotes the unit square (2D case) or cube (3D case); the coefficients p, q
and r are positive and bounded while t is nonnegative and bounded. The differential
operator is discretized by using the five-point (2D case) or seven-point (3D case) finite
difference approximation (point scheme box integration [30, 41]). The mesh points are
ordered according to the lexicographic ordering. One then results in a large system
of linear equations of the form
where A is a block tridiagonal or pentadiagonal diagonally dominant Stieltjes matrix,
b is a vector that depends on both the rhs f of (1.1) and the boundary conditions,
while u is the vector of unknowns. Combined with some appropriate preconditioning
matrix, the conjugate gradient method is (one of) the most widely used method(s) (see
e.g. [1, 2, 6, 11, 19, 20, 35]) for solving system (1.2). Relaxed incomplete factorizations
(RILU) are powerful preconditioning techniques that interpolate, through a relaxation
parameter, between the popular incomplete LU factorization ILU and its modified
variant MILU that preserves rowsums of A [4, 13]. As opposed to ILU and MILU,
the two main advantages of RILU are following
y Research supported by the Commission of the European Communities HCM Contract No. ERB-
CHBG-CT93-0420, at Utrecht University, Mathematical Institute, The Netherlands. Current address
de Bruxelles, Service des Milieux Continus (CP 194/5), 50, avenue F.D.
Roosevelt, B-1050 Brussels, Belgium. magolu@ulb.ac.be
z Universit'e Libre de Bruxelles, Service de M'etrologie Nucl'eaire (CP 165), 50 av. F.D. Roosevelt,
B-1050 Brussels, Belgium. ynotay@ulb.ac.be. Research supported by the " Fonds National de la
Recherche Scientifique", Belgium.
M.M. MAGOLU AND Y. NOTAY
1. it does not suffer a lot from existence problem [13, 17, 18];
2. it is robust with respect to discontinuities and anisotropy [4, 34].
Two major inconveniences are that the "optimal" value of the relaxation parameter
strongly varies from a problem to another and the behavior could be very
sensitive to variations of ! around the observed "! opt " [12, 39]. In [34], a new variant
of RILU has been proposed. There, the relaxation parameter is variable and dynamically
computed during the incomplete factorization phase. Like its precursor RILU,
it is robust with respect to both existence and performance. In addition, its performance
does not critically depends on involved parameter. We intend to propose a
block version of this dynamically relaxed method (DRBILU) and to compare it with
classical pointwise and blockwise methods as well as with some existing "dynamic"
variants. We stress that, even for three-dimensional problems, we consider a linewise
partitioning of the unknowns : each block corresponds to a set of gridpoints along a
line parallel to the x-axis in the physical domain [28].
Our study is outlined as follows. Needed general terminology and notation are
gathered in Section 2. In Section 3, we first review some variants of block incomplete
factorizations in Subsection 3.1 and 3.2. We next establish, in Subsection 3.3,
some theoretical results that motivate the introduction of dynamically relaxed block
preconditioners. Comparative numerical experiments are reported and discussed in
Section 4.
2. General terminology and notation.
2.1. Order relation. The order relation between real matrices and vectors is
the usual componentwise order : if
A is called nonnegative (positive) if A - 0
2.2. Stieltjes matrices. A real square matrix A is called a Stieltjes matrix (or
equivalently, a symmetric M-matrix) if it is symmetric positive definite and none of
its offdiagonal entries is positive [10].
2.3. Normalized point LU-factorization. Given a Stieltjes matrix S, by its
normalized point LU factorization we understand the factorization
s where
P s is pointwise diagonal and L s is pointwise lower triangular such that diag(L s
2.4. Miscellaneous symbols. We describe below some symbols that are used
in our study. A denotes a given square matrix of order n.
the transpose of A
the ith smallest eigenvalue of A
(A), the smallest eigenvalue of A
(A), the largest eigenvalue of A
the pointwise diagonal matrix whose diagonal
entries coincide with those of A
the pointwise tridiagonal matrix whose main
three diagonals coincide with those of A
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 3
diag the block diagonal matrix whose block
diagonal entries coincide with those of A
offdiag
e : the vector whose all components are equal to 1
Throughout this work, the term block will refer to the linewise partitioning mentioned
in the introduction.
3. Blockwise incomplete factorizations. For simplicity, throughout this work,
we consider only blockwise incomplete factorizations with no fill-in allowed outside
the main block diagonal part of A. Let n x , n y and n z denote the number of unknowns
in respectively the x-, the y- and the z-direction (if any), then the order of the matrix
A is According to our assumptions, A is either block tridiagonal (2D
case), i.e.
or block pentadiagonal (3D case), i.e.,
A may be split as
into its block diagonal part D and strictly block lower (upper) triangular part \GammaL
(\GammaL t ). The matrix
where P is the block diagonal matrix computed according to Algorithm 3.1, is referred
to as the block incomplete LU-factorization (BILU) of A [16]. P
stands for the normalized point LU factorization of P i;i . Other choices for g
are
discussed in [8, 16]. Algorithms that handle more general matrices may be found in
[8, 25].
Given that all the blocks A i;i\Gamma1 and A i;i\Gamman y are diagonal, each P i;i is a tridiagonal
Stieltjes matrix. It is well known that, if T is a pointwise tridiagonal Stieltjes
may be cheaply computed from the normalized point LU
factorization of
T . This follows from the fact that T \Gamma1 may be rewritten
as
The lower part of e
may then be computed by simple identification of the relevant
entries in the above identity. This is performed in Algorithm 3.2 where, for simplicity,
we set
To improve the performance of Algorithm 3.1, various variants have been proposed
in the literature. Most of them differ from the basic method only in the way of
computing diag(P i;i ) which is modified so as to satisfy rowsum relations and/or to keep
control of the extreme eigenvalues of the preconditioned matrix B \Gamma1 A [4, 21, 22, 26].
Reviewing all existing block variants is beyond the scope of our study. We shall
confine ourselves to relaxed and dynamically modified methods.
4 M.M. MAGOLU AND Y. NOTAY
Algorithm 3.1 (BILU)
Compute, for
if (3D) and ny
i\Gamman y ;i\Gamman y
A i\Gamman y ;i
(2)
LU-factorization
Algorithm 3.2 ( e
qn;n
normalized point LU-factorization of T
3.1. Relaxed block incomplete factorizations. Now, fill-ins that are neglected
in Algorithm 3.1 are accumulated and added to the current main diagonal
after multiplication by a relaxation parameter !
3.3. It amounts to impose the following rowsum relation [22]
where M stands for the pointwise diagonal matrix defined as in Algorithm 3.3 with
In other words, M satisfies
diag block
In the 2D case, (3.4) simplifies to
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 5
Algorithm 3.3 (RBILU)
Compute, for
with
if (3D) and ny
i\Gamman y ;i\Gamman y
A i\Gamman y
with
A i;i\Gamman y
s i\Gamman
i\Gamman y ;i\Gamman y A i\Gamman y ;i e i
(2)
if (3D) and
A i;i+1 e i+1 +A i;i+ny e i+ny
relaxed block incomplete LU-factorization
diagonal (modification) matrix
As is well known, multiplication by P \Gamma1 (or by P
in Algorithm 3.3) has to
be implemented as solution of a (small) linear system. With recovers BILU
while with gets the standard modified variant (MBILU) whose conditioning
properties are theoretically investigated in, a.o., [9, 27, 32, 23, 5, 25, 28]. In the case
under consideration, A is a (nonsingular) irreducible Stieltjes matrix, Algorithm 3.3
cannot breakdown; it gives rise to a (nonsingular) diagonally dominant Stieltjes matrix
later purposes, note that the existence analysis discussed in [26, 32]
covers any incomplete factorization such that
Be
Like their pointwise counterparts, both standard methods
suffer mainly from robustness problems [17, 22, 25, 26, 34, 40]. Optimal performances
are achieved with 0- 1. The trouble is that ! opt strongly depends on the problem
Most severe is the fact that performances could be highly sensitive to the
variation of ! around ! opt which is very hard to estimate [4, 12, 39]. In the case of
6 M.M. MAGOLU AND Y. NOTAY
uniform grid of mesh size h in all directions, it is advocated in [4] to use
problems [4]), in which case one has (see Subsection 3.3) that
In the light of the theory extensively developed in [9, 22, 23, 28, 32, 34], one should
take
with
ae n y in 2D case
case
in order to handle any grid. In [39], it has been suggested to try In [34] a
pointwise dynamic version of RILU, termed DRILU, that improves the performance
stability with respect to the parameter involved, has been introduced. Before discussing
the blockwise version of DRILU, we would like to say a few words about
dynamically modified block methods.
3.2. Dynamically modified block incomplete factorizations. The standard
modified method efficient only in "nice" situations; e.g., in the case
of fixed mesh size, Dirichlet boundary conditions and monotonous variation for the
PDE coefficients. It is now well-established that the performance strongly depends on
the ordering strategy, the variations of both the PDE coefficients and the mesh size,
and the boundary conditions [25, 28]. With dynamically modified methods, the goal
is to be in more general situations as efficient as, or even more efficient than, classical
modified method in nice circumstances, without changing the numbering of the
unknowns. To this end, small perturbations are dynamically added to the diagonal
entries of P i;i (initially computed with imposed constraints are
stand for a O(1) positive parameter independent of n x , n y and n z . If
one applies Algorithm 3.4 then the following rowsum relation holds
where, at each grid point j, the perturbation - j is defined as in Algorithm 3.4(3).
There holds Moreover, one has that [22, 26]
In [22, 26], where the 2D case is discussed, it is proposed to take
we extend to the 3D case as in agreement with the theoretical argument
in [9, 22, 23, 28, 32]. The analysis of the 2D case shows that the performances do not
strongly depend on the variation of
and that "i opt - 1
range of PDEs [22, 26].
Algorithm 3.4 outperforms by far both basic block methods (Algorithm 3.3 with
long as there is isotropy (including strong discontinuities) or moderate
anisotropy in the PDE coefficients (see [22]). In the case of strong anisotropy, it
gives rise to several degenerate (say, of O(fflh 2 ) with isolated smallest
eigenvalues, which slows down the convergence of PCG process [33, 37, 38, 26]. This
occurs e.g. in the case of PDEs that involve both isotropy and strong anisotropy (see
e.g. [24, 26, 34]). In such a situation, in the 2D case, it is better to use Algorithm 3.5
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 7
Algorithm 3.4 (DMBILU)
Compute, for
with
if (3D) and ny
i\Gamman y ;i\Gamman y
A i\Gamman y
with
+A i;i\Gamman y
s i\Gamman
i\Gamman y ;i\Gamman y A i\Gamman y ;i e i
e
if (3D) and
e
A
e
dynamically modified block incomplete LU-factorization
diagonal (modification) matrix
that cancels the perturbation - j at unsafe nodes [24, 26]. Observe that Algorithm 3.5
reduces to Algorithm 3.4 in the absence of "strong" anisotropy. Now there holds
where, at each grid point j, the perturbation - j is equal to - j whenever - j is defined,
and 0 otherwise [24, 26]. It is worth noting that unsafe nodes are nodes where the
coefficients p and q are strongly anisotropic; the factor 10 in Algorithm 3.5 gives
a "measure" of the amount of anisotropy beyond which the perturbations should be
discarded [26, Subsection 4.4]. Not any generalization of Algorithm 3.5 has been
proposed so far to handle 3D problems. This is not a trivial task; one has to take into
account the PDE coefficient r too. Contrary to the 2D case [24, 26], it is not easy to
8 M.M. MAGOLU AND Y. NOTAY
establish in which of the many different possibilities, with respect to anisotropy, one
should drop the perturbations without dramatically increasing the largest eigenvalues
of the preconditioned matrix. An alternative solution is investigated in the next
section, which is the main contribution of this work.
Algorithm 3.5 (DMBILU   ) (2D case)
Compute, for
as in DMBILU
then
endif
as in DMBILU
dynamically modified block incomplete
LU-factorization with dropping test
3.3. Dynamically relaxed block incomplete factorizations. Relaxed methods
have been observed to successfully handle strongly anisotropic problems, whenever
the relaxation parameter is properly chosen (see, e.g., Section 4). Unfortunately, no
general theory has been provided to date to estimate ! opt . The pointwise variant of
standard RILU, introduced in [34] and termed DRILU, dynamically computes variable
relaxation parameters; it combines the robustness, with respect to the variation of the
PDE coefficients, of optimized RILU with the relative insensitivity, to the variation of
the parameters involved, of dynamically perturbed methods. The blockwise version
of DRILU that we are going to present in this section is based on a generalization
of [34, Theorem 5.1] which motivated the introduction of DRILU. We first give an
auxiliary result.
Lemma 3.1. [34, Lemma 5.1] Let F - 0 be a pointwise strictly upper triangular
matrix and, Q 0 and Q 1 stand for nonnegative pointwise diagonal matrices. If C is a
matrix and x a positive vector such that
then C is nonnegative definite.
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 9
Theorem 3.2. Let stand for a diagonally dominant Stieltjes
matrix, such that while L is strictly block lower triangular. Let P
denote a block diagonal diagonally dominant Stieltjes matrix.
p be the normalized point LU factorization of P . Assume
further that
diag block
tridiag
M is a pointwise diagonal matrix such that
f
diag block
e
with
for some pointwise diagonal matrix
such that W - I.
If, for all 1
where
then
Proof. First, given that
to check that, (3.12) and (3.13) imply that
Be
diag block
e
offdiag block
e
diag block
offtridiag
Next, let be the pointwise diagonal matrix whose diagonal entries are
given by
One has 0 - \Theta - \Gamma1 I. Set
diag block
tridiag
denotes the pointwise diagonal matrix such that B 1 Taking (3.12)
into account, it is an easy matter to show that offdiag(B 1 Hence it follows that
is a Stieltjes matrix and therefore nonnegative definite [10]. Now,
p is
a Stieltjes matrix such that P e - 0, whence (see [10]) L \Gamma1
be the pointwise diagonal matrix defined by
ae
one has by (3.16) that XP2
denote the pointwise diagonal matrices such that, respectively,
By application of Lemma 3.1, successively to B 2 and B 3 , one
easily shows that both matrices are nonnegative definite. On the other hand, since
readily checks that
diag block
offtridiag
offdiag block
where \Delta 4 is a pointwise diagonal matrix whose explicit form does not matter. By
(3.18) one has that
diag block
offtridiag
e
offdiag block
Now, by definition of - (see (3.17)), one easily deduces that 1+! i
i, so that either -
This implies that I \Gamma W - 2\Theta. Therefore the right hand side of (3.19) is nonnegative
definite. The conclusion readily follows.
Corollary 3.3. If B corresponds to some relaxed block incomplete LU factorization
of A with
with
Proof. Apply Theorem 3.2 with M defined as in Eq. (3.4) and
all
Note that, since L \Gamma1
imply for MBILU factorizations
with
(3.
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 11
(3.22) is nothing but the upper bound on the basis of which both DMBILU and
DMBILU   have been elaborated, by imposing ~
- i' [22, 24, 26]. An alternative way
to achieve the latter imposed upper bound consists in computing P so as to satisfy
the assumptions of Theorem 3.2 with
for all
In view of (3.12) and (3.13), this could be achieved by
1. computing the entries of P as in Algorithm 3.1 (block by block);
2. subtracting, from the pointwise diagonal entries of P , the quantity ( f
which is defined by (3.13) with equality symbol, where ! i is defined by (3.24).
The corresponding block preconditioner, which we call dynamically relaxed block incomplete
LU-factorization (DRBILU), is described in Algorithm 3.6. Note that, if
p is a pointwise tridiagonal Stieltjes matrix, then so is
.
Therefore, e
may be computed by means of Algorithm 3.2 with
Unlike RBILU, the relaxation parameter is now dynamically
modulated in function of - i 's, in a similar way as perturbations are added in
DMBILU. This leads us to expect a similar stability with respect to the choice of the
parameter i.
For all the block preconditioners discussed so far, the parameters involved may be
chosen so as to achieve the same upper bound, say i', for the the largest eigenvalues of
the preconditioned matrix B \Gamma1 A. In the context of the PCG method, the convergence
behavior also depends on the distribution of the smallest eigenvalues. The following
estimate, that relates - has been obtained
and successfully tested in [7, 26, 31]
where X denotes the pointwise diagonal matrix such that Be
which satisfies \Gamma1 - ffl h;i - 1, depends weakly upon both the mesh size parameter h
and i. ffl rise in general to very accurate estimates in the case of i - n,
for both pointwise [7, 31] and blockwise preconditioners [26]. One has (e;
that the order of
magnitude of the smallest eigenvalues - mostly depends on the sum (e; X e) of
all perturbations. As far as DMBILU and DMBILU   are concerned one has (e; X
O(1) [22, 26], whence it follows that the smallest eigenvalues of B \Gamma1 A do not depend
on h, i.e., are O(1).
As regards RBILU, with
i' , one has from (3.3) that
i' M , which
shows that all perturbations but the first block nodes ones are O(h). The smallest
eigenvalues are then clearly O(h). Therefore, -(B \Gamma1 asymptotically O(h \Gamma2 ).
Nevertheless, with a "properly chosen" value for !, RBILU performs very well in
practice because of the nice distribution of (interior) eigenvalues [4].
Now, for DRBILU where
perturbations (i.e. occur only at
selected nodes according to (3.24). Moreover, for all
i' . It is then an easy matter to show that, for the same target
upper bound, the perturbations introduced in DRBILU are not larger than the ones
added in RBILU. Therefore, on the basis of (3.25) one may expect DRBILU to be at
least as robust as RBILU. It is of worth noting that for both RBILU and DRBILU,
Algorithm 3.6 (DRBILU)
Compute, for
with
e
if (3D) and ny
i\Gamman y ;i\Gamman y
A i\Gamman y
with
e
+A i;i\Gamman y
s i\Gamman
i\Gamman y ;i\Gamman y A i\Gamman y ;i e i\Gamman y
(2)
e
if (3D) and
e
dynamically relaxed block incomplete LU-factorization
diagonal (relaxation) matrix
the perturbations are in direct proportion to the neglected fill-ins which are (very)
small; this occurs in particular in 2D problems with strong anisotropy (p - q or
Subsection 4.4] (see also (3.28)).
For comparison purposes, let us mention that as far as DMBILU is concerned,
one has by [22, Lemma 4.4], with
O(h), the following sharp upper bound
on the perturbations -
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 13
This bound could be large in the case of strong anisotropy. For instance, in 2D
case, when p - q (in (1.1)) one could have that - In
such a situation, the smallest degenerate eigenvalues of D \Gamma1 A are reproduced, up
to some multiplicative constants, by B \Gamma1 A (see (3.25)). If the number - of such
degenerate eigenvalues is not very small (i.e., if
n), this results in a very
slow convergence for the PCG process [33]. In the case of DRBILU, since (3.16) is
equivalent to - i
, it is straightforward to establish that
with
Therefore,
1-j-n
diag block
tridiag
e
whence it follows that the perturbations involved in DRBILU can never be much
larger than the ones involved in DMBILU. It is obvious from all the considerations
above that, with DRBILU, one tries to combine the advantages of both RBILU and
DMBILU.
Observe finally that, in the case of 2D problems, even if the perturbations introduced
by DRBILU are very small, their sum (e; X
Me) could be larger
than the corresponding sum for DMBILU   , in particular when the number j of nodes
where the perturbations are dropped is large enough, for instance O(n). In order
to perform a fair comparison, we give in Algorithm 3.7 a version of DRBILU, termed
DRBILU   , which uses the same dropping test as in DMBILU   .
4. Numerical experiments. The PCG method is run with the zero vector as
starting approximate solution and the residual error reduction
as convergence criterion. The computations are performed in double precision FORTRAN
on a Sun 514MP sparc workstation. For comparison purposes, the precondi-
tionings include :
1. RBILU (Algorithm 3.3). Four values of the parameter ! have been tested :
respectively, the unmodified and (unper-
turbed) modified standard block methods;
defined as in Eq. (3.8). For small
and moderate problem sizes, this includes ! - 0:95 that has been suggested
in [39], while for large problems, this includes ! - 0:99 which is the optimal
observed to the nearest 0:01 for minimizing the number of iterations.
14 M.M. MAGOLU AND Y. NOTAY
Algorithm 3.7 (DRBILU   ) (2D case)
Compute, for
as in DRBILU
then
else
endif
as in DRBILU
dynamically relaxed block incomplete LU-factorization
with dropping test
2. DMBILU and DMBILU   (Algorithms 3.4 and 3.5). We have used 1according to the recommendations made in [26].
3. DRBILU and DRBILU   (Algorithms 3.6 and 3.7). We report the results
4 which we anticipate to be near optimal as in DMBILU and
DMBILU   .
4. ILU and MILU. The popular pointwise unmodified and modified incomplete
LU-factorizations (see, e.g., [6]).
5. DRILU ([34]). The pointwise dynamically relaxed incomplete LU-factorization
method. As recommended in [34], the parameter involved is chosen so as to
target the upper bound n 1
d where denotes the space dimension.
In the first two problems, DMBILU   and DRBILU   are not considered because
they coincide with, respectively, DMBILU and DRBILU. Problem 3 is essentially
Stone's problem [36, 26]. The next three problems are some 3D extensions of the first
three ones. The last example is intended for comparing the behavior of block and
pointwise methods in the case of elongated grids and non uniform mesh size, which
arise in 3D simulation problems.
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 15
Problem 1. (2D)
ffl The rhs of the linear system to solve is chosen such that the function u 0 (x;
generates the solution on the grid.
Problem 2. (2D, from [40])
ae 100 in (1=4; 3=4) \Theta (1=4; 3=4)
elsewhere
ae 100 in (1=4; 3=4) \Theta (1=4; 3=4)
elsewhere
Problem 3. (2D, essentially from [36])
ffl The coefficients p, q and t are specified in Fig. 1
ffl The rhs of the linear system to solve is chosen such that the function u 0 (x;
generates the solution on the grid.
-x
y
region p q
d
d
t0Fig. 1. Problem 3. Configuration and specification of the PDE coefficients; d stands for a positive
parameter.
Problem 4. (3D)
@\Omega
M.M. MAGOLU AND Y. NOTAY
Problem 5. (3D)
ae 100 in (1=4; 3=4) \Theta (1=4; 3=4) \Theta (0; 1)
elsewhere
ae 100 in (1=4; 3=4) \Theta (1=4; 3=4) \Theta (0; 1)
elsewhere
Problem 6. (3D)
ffl The coefficients p, q, r, t and f depend only on (x,y) as specified in Fig. 2
@\Omega -x
y
region p q
d
d
r1d
Fig. 2. Problem 6. Configuration and specification of the PDE coefficients; all the regions extend
from stands for a positive parameter.
Problem 7. (3D)
ae
1000 in (0; 1=8) \Theta (0; 1=8) \Theta (0; 1=8)
elsewhere
ae 1 in (0; 1=8) \Theta (0; 1=8) \Theta (0; 1=8)
elsewhere
For all problems but the last one, we have used a uniform mesh size h in each
direction. In the case of Problem 7, we have used non uniform rectangular grids
obtained by setting h
ny and h z = 4h0 (z)
nz , where the function h 0
is defined by
0:25 if 0:25 - t - 0:5
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 17
We give in Tables 1-14 the extremal eigenvalues and/or the spectral condition
numbers as well as the exponent - from the assumed asymptotic relationship -(B \Gamma1
constant. - is estimated from the largest two problems
data (say, h \Gamma1 =96 and h \Gamma1 =192, or 40 and 80 in 3D cases). Whenever - min (B \Gamma1
is strongly isolated from the rest of the spectrum, we also include both the second
smallest eigenvalue and the effective spectral condition number
which accounts for the superlinear convergence of the PCG method (see e.g. [33,
37, 38, 40]). As far as pointwise preconditioners (ILU, MILU and DRBILU) are
concerned, whose conditioning analysis is not investigated here, we have computed
only the number of iterations to achieve the prescribed accuracy, in order to save
space. Problems 3 and 6, with d large enough, are examples of PDE with degenerate
smallest eigenvalues. We report in Tables 5 and 11, for
case) and h case), the numerically computed smallest and largest four
eigenvalues associated to each block preconditioner involved. The pointwise Jacobi
(or diagonal) preconditioner whose smallest eigenvalues are connected to those of the
other preconditioners through Eq. (3.25), is also considered. In Tables 6, 13 and 14,
the number of PCG iterations to reach the target accuracy are collected for each
problem and for h
of the mesh sizes along the three directions in the case of Problem 7. From all the
tables, the following observations can be made.
1. As expected, RBILU (with smallest eigenvalues of O(h)
and largest eigenvalues of approximately O(h \Gamma1 ). Nevertheless, the O(h \Gamma2 )
behavior of -(B \Gamma1 compensated by the nice distribution of interior eigen-values
[4], which explains the relative good performance of the preconditioner
(see

Tables

6, 13 and 14, RBILU with
2. For both DMBILU and DRBILU as well as their "  versions", the smallest
eigenvalues are in general O(1) while the largest ones are O(h \Gamma1 ), so that the
(effective) spectral condition numbers are O(h \Gamma1 ). Observe however that, in
the case of strong anisotropy (Problems 3 and 6 with smallest
eigenvalues associated to DRBILU are O(h), whence it follows that the
(effective) spectral condition numbers are O(h \Gamma2 ); as in RBILU, the good
behavior of DRBILU is due to the nice distribution of interior eigenvalues.
3. In the case of isotropic problems, DMBILU and DRBILU give better results
than classical methods (RBILU with does "optimized" RBILU.
All winning methods perform quite similarly, even if there is strong jump
in the PDE coefficients (Problems 2 and 5), in which case the small first
eigenvalue is the only one that is strongly isolated from the others.
4. In the presence of strong anisotropy (Problems 3 and 6 with d=10 3 ), RBILU
essentially reproduce the disastrous distribution of
the smallest eigenvalues of D \Gamma1 A, which slows down the convergence of PCG
process [33]. It should be noted that for RBILU with
eigenvalues are more clustered than for DMBILU as one would expect by
comparing the largest eigenvalues. DRBILU, DMBILU   and DRBILU   successfully
break with the dependence upon D \Gamma1 A, which reflects in the number
of PCG iterations. As predicted, DRBILU   is a little bit more efficient than
DRBILU. For RBILU with the rate of convergence mostly depends on
M.M. MAGOLU AND Y. NOTAY
the distribution of the largest eigenvalues, the smallest ones being known to
cluster around 1.
5. In accordance with previous works, [16] (2D), [3, 28] (3D), blockwise (linewise)
methods turn out to be more efficient than pointwise counterparts. In the
case of 2D problems, the reduction in the number of iterations, from point
methods to block methods, is at least about 50%, while it is around 30%
in the 3D cases. The gain is even more spectacular in the case of strongly
anisotropic problems (see Table 6, Problem 3b and Table 13, Problem 6b).
As regards their computational complexity, let us mention that each PCG
iteration with blockwise preconditioners needs two more flops per point than
for pointwise preconditioners, which is rather small as compared to the total
number of flops per PCG iteration [3].
6. The performances of blockwise methods in general, and in particular, DM-
BILU and DRBILU, are (almost) insensitive to the variation of the number
of gridpoints along the x-direction, say, the direction which determines the
blocks (see Tables 12 and 14). This is in quite agreement which the analysis
performed in [28].
7. DRBILU is the only preconditioner which is always among the best three
ones, whatever the problem tested (see Tables 6, 13 and 14). Whenever
DRBILU is not absolutely the winner, it is not far from the latter (we have
also observed that the variation of the parameter i around 1
4 does not have
a significant effect on the behavior of the preconditioners).
From our discussion together with the analysis made in Section 3, we conclude
that the most promising method is DRBILU (with It perfoms quite well
in a wide range of situations. Its main merit is that, contrary to DMBILU and
DMBILU   , it does not (strongly) depend on a dropping test that has been set up on
the basis of experiments performed only on five-diagonal 2D problems [24, 26]. The
dropping test involved has not yet been extended to 3D PDEs, while the performances
of DRBILU are quite satisfactory for both two- and three-dimensional problems. Even
though all our computations were performed on well structured grids for which natural
blockwise partitionings are available, the block methods that we have investigated
could be applied to finite element unstructured grids. The problem of defining blocks
in such irregular grids has been solved recently in [14, 15]. Combining the techniques
developed in the latter papers with the preconditoners discussed here would result in
robust preconditioners to tackle real life engineering problems. This awaits further
investigation.
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 19

Table
Problem 1. Extremal eigenvalues (-min and -max) and/or spectral condition number (-) of B
exponent - corresponding to the (estimated) asymptotic relationship -=Ch \Gamma- , C denoting a constant.
RBILU
48 16.2 4.02 0.58 2.39 4.09 0.39 1.96 5.05
DMBILU DRBILU
48 0.65 2.52 3.88 0.77 2.80 3.64 0.98 3.78 3.84

Table
Problem 2. Extremal eigenvalues (-min , -2 and -max) and/or (effective) spectral condition number
corresponding to the (estimated) asymptotic relationship
denoting a constant.
RBILU
48 2300 13.7 51.3 70E-4 0.61 2.97 422 4.86
48 38E-4 0.44 2.28 608 5.21 54E-4 0.56 2.94 549 5.27
DRBILU
48 76E-4 0.66 3.29 435 4.98 20E-3 0.85 5.42 271 6.39
M.M. MAGOLU AND Y. NOTAY

Table
Problem 3 with Extremal eigenvalues (-min , -2 and -max) and/or (effective) spectral
condition number (- (2) ) - of B \Gamma1 A; exponent - corresponding to the (estimated) asymptotic relationship
-=Ch \Gamma- , C denoting a constant.
RBILU
48 7298 38.4 65.8 25E-4 0.42 3.54 8.52
192 116422 614. 309. 57E-5 0.11 6.85 64.0
48 13E-4 0.24 2.76 11.5 50E-5 0.13 3.22 24.8
DRBILU
48 19E-4 0.34 3.63 10.6 44E-4 0.67 5.98 8.95
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 21

Table
Problem 3 with Extremal eigenvalues (-min , -2 and -max) and/or (effective) spectral
condition number (- (2) ) - of B \Gamma1 A; exponent - corresponding to the (estimated) asymptotic relationship
-=Ch \Gamma- , C denoting a constant.
RBILU
48 7302 36.2 75.0 26E-4 0.43 3.14 7.22
48 14E-4 0.28 2.45 8.88 11E-6 37E-4 3.84 1029
DRBILU
48 19E-4 0.34 4.20 12.4 43E-4 0.64 6.49 10.2
22 M.M. MAGOLU AND Y. NOTAY

Table
Problem 3 with Distribution of extremal eigenvalues of B \Gamma1 A for different
preconditioners.
Preconditioning smallest eigenvalues largest eigenvalues
1. 1. 1. 1. 159. 204. 248. 368.
Point Jacobi 3E-9 59E-8 17E-7 46E-7 2. 2. 2. 2.

Table
Number of PCG iterations to achieve kr (i)
Problem 1 Problem 2 Problem 3a Problem 3b
ILU
DYNAMICALLY RELAXED BLOCK PRECONDITIONERS 23

Table
Problem 4. Extremal eigenvalues (-min and -max) and/or spectral condition number (-) of B
exponent - corresponding to the (estimated) asymptotic relationship -=Ch \Gamma- , C denoting a constant.
RBILU
DMBILU DRBILU

Table
Problem 5. Extremal eigenvalues (-min , -2 and -max) and/or (effective) spectral condition number
corresponding to the (estimated) asymptotic relationship
denoting a constant.
RBILU
DRBILU
M.M. MAGOLU AND Y. NOTAY

Table
Problem 6 with Extremal eigenvalues (-min , -2 and -max) and/or (effective) spectral
condition number (- (2) ) - of B \Gamma1 A; exponent - corresponding to the (estimated) asymptotic relationship
-=Ch \Gamma- , C denoting a constant.
RBILU
DRBILU

Table
Problem 6 with Extremal eigenvalues (-min , -2 and -max) and/or (effective) spectral
condition number (- (2) ) - of B \Gamma1 A; exponent - corresponding to the (estimated) asymptotic relationship
-=Ch \Gamma- , C denoting a constant.
RBILU
DRBILU

Table
Problem 6 with Distribution of extremal eigenvalues of B \Gamma1 A for different
preconditioners.
Preconditioning smallest eigenvalues largest eigenvalues
1. 1. 1. 1. 444. 511. 807. 3129
Point Jacobi 17E-9 39E-7 57 E-7 15E-6 2. 2. 2. 2.
26 M.M. MAGOLU AND Y. NOTAY

Table
Problem 7. Extremal eigenvalues (-min and -max) and/or spectral condition number (-) of B \Gamma1 A.
RBILU
grid -min -max -min -max -
160 \Theta 80 \Theta 40 374. 30.2 0.12 10.6 91.7 62E-3 7.5 121.
DMBILU DRBILU
grid -min -max -min -max -min -max -
160 \Theta 80 \Theta 40 0.29 10.6 36.7 0.33 11.0 32.4 0.63 15.3 24.1

Table
Number of PCG iterations to achieve kr (i) Problems 4, 5, and 6 (6a:
Problem 4 Problem 5 Problem 6a Problem 6b
ILU
DRILU

Table
Number of PCG iterations to achieve kr (i) k2 =kr (0) k2 -10 \Gamma7 for Problem 7. Grids: (a)=40 \Theta 40 \Theta
(g)=160 \Theta 80 \Theta 40 , (h)=80 \Theta 160 \Theta 40 .
grid (a) (b) (c) (d) (e) (f) (g) (h)
43 43 61
28 42 28 28 28
ILU 76 158 114 128 175 138 205 205

Acknowledgments

. Part of this work was done while the first author was holding
a postdoctoral position at the Mathematical Institute of Utrecht University. He
thanks Henk van der Vorst for his warm hospitality, and for suggesting to include the
three-dimensional case in this study. Thanks are also due to the referees for their
constructive comments.



--R

Cambridge University Press
Finite Element Solution of Boundary Value Problems.
Vectorizable preconditioners for elliptic difference equations in three space dimensions
On the eigenvalue distribution of class of preconditioning methods
On eigenvalue estimates for block incomplete factorization methods
Templates for the Solution of Linear Systems
Modified incomplete factorization strategies
On sparse block factorization iterative methods
Existence and conditioning properties of sparse approximate block factorizations
Nonnegative Matrices in the Mathematical Sciences
A survey of preconditioned iterative methods
Fourier analysis of relaxed incomplete factorization preconditioners
Approximate and incomplete factorizations
An object-oriented framework for block preconditioning
BPKIT block preconditioning tool kit
Block preconditioning for the conjugate gradient method
Beware of unperturbed modified incomplete factorizations
Relaxed and stabilized incomplete factorizations for nonself-adjoint linear sys- tems
Matrix Computations
Closer to the solution: iterative linear solvers
Compensative block incomplete factorizations
Modified block-approximate factorization strategies
Analytical bounds for block approximate factorization methods
Empirically modified block incomplete factorizations
Ordering strategies for modified block incomplete factorizations
Taking advantage of the potentialities of dynamically modified block incomplete factorizations
On the conditioning analysis of block approximate factorization methods
Theoretical comparison of pointwise
Efficient planewise like preconditioners to cope with 3D prob- lems
Computational Methods in Engineering and Science

Conditioning analysis of modified block incomplete factorizations
On the convergence rate of the conjugate gradients in the presence of rounding errors

Iterative Methods for Sparse Linear Systems
Iterative solution of implicit approximation of multidimensional partial differential equations
The convergence behaviour of conjugate gradients and ritz values in various circumstances
The rate of convergence of conjugate gradients
ICCG and related methods for 3D problems on vector computers
The convergence behaviour of preconditioned CG and CG-S
Iterative Solution of Elliptic Systems and Application to the Neutron Diffusion Equations of Reactor Physics
--TR
large sparse linear systems;preconditioned conjugate gradient;incomplete factorizations;diagonal relaxation;discretized partial differential equations