Schur Complement Systems in the Mixed-Hybrid Finite Element Approximation of the Potential Fluid Flow Problem.
The mixed-hybrid finite element discretization of Darcy's law and continuity equation describing the potential fluid flow problem in porous media leads to a symmetric indefinite linear system for the pressure and velocity vector components. As a method of solution the reduction to three Schur complement systems based on successive block elimination is considered. The first and second Schur complement matrices are formed eliminating the velocity  and pressure variables, respectively, and the third Schur complement matrix is obtained by elimination of a part of Lagrange multipliers that come from the hybridization of a mixed method. The structural properties of these consecutive Schur complement matrices in terms of the discretization parameters are studied in detail. Based on these results the computational complexity of a direct solution method is estimated and compared to the computational cost of the iterative conjugate gradient  method applied to Schur complement systems. It is shown that due to special block structure the spectral properties of successive Schur complement matrices do not deteriorate and the approach based on the block elimination and subsequent iterative solution is well justified. Theoretical results are illustrated by numerical experiments.
Introduction
.
Let
be a bounded domain in R 3 with a Lipschitz continuous
boundary @
The potential
uid
ow in saturated porous media can be described by
the velocity u using Darcy's law and by the continuity equation for incompressible
ow
where p is the piezometric potential (
uid pressure), A is a symmetric and uniformly
positive denite second rank tensor of the hydraulic resistance of medium with
for all represents the density of potential sources in the
medium. The boundary conditions are given by
@
@
where
@
@
@
N are such that
@
@
@
is the outward
normal vector dened (almost everywhere) on the boundary @
Assume that the
domain
is a polyhedron and it is divided into a collection of
subdomains such that every subdomain is a trilateral prism with vertical faces and general
nonparallel bases (see, e.g., [11], [14] or [15]). We will denote the discretization of the
domain
by E h and assume an uniform regular mesh with the discretization parameter
h. Denote also the collection of all faces of elements which are not adjacent
This work was supported by the Grant Agency of the Czech Republic under grant 201/98/P108 and
by the grant AS CR A2030706. Revised version October 1999.
y Seminar for Applied Mathematics, Swiss Federal Institute of Technology (ETH) Zurich, ETH-
Zentrum, CH-8092 Zurich, Switzerland. (miro@sam.math.ethz.ch)
z Institute of Computer Science, Academy of Sciences of the Czech Republic, Pod vodarenskou vez
2, 182 07 Prague 8, Czech Republic. (maryska@uivt.cas.cz, miro@uivt.cas.cz, tuma@uivt.cas.cz)
to the boundary
@
D by
@
D and introduce the set of interior faces
@
We consider the following low order nite element approximation. Let
be the space spanned by the linearly independent basis functions v e
dened on the element e 2 E h in the form
and such that they are orthonormal with respect to the set of functionals
Z
5:
Here f e
k denotes the k-th face of the element e
k;3 ) is the
outward normal vector with respect to the face f e
k . The velocity function u will be
then approximated by vector functions linear on every element e 2 E h from the Raviart-
Thomas space
denotes the restriction of a function v h onto the element e 2 E h . Further
denote the space of constant functions on each element e 2 E h by M 0 (e) and denote the
space of constant functions on each face f 2 h by M 0 (f ). The piezometric potential p
will be approximated by the space which consists of elementwise constant functions
where  h j e is the restriction of a function  h onto element e . The Lagrange multipliers
coming from the hybridization of a method will be approximated by the space of
all functions constant on every face from h
Here  h j f denotes the restriction of a function  h onto the face f 2 h . Analogously we
introduce the spaces M 0(@
N ) as the spaces of functions constant on
every face from [ e2E h
@e \
@
D and h \
@
respectively. The detailed description of
the spaces that we use can be found in [14] (see also [11] or [15]).
The Raviart-Thomas approximation of the mixed-hybrid formulation for the problem
(1.1) and (1.2) reads as follows (see [4]):
Find
@e\@
(1.
@e\@
where p D;h and uN;h are approximations to the functions p D and uN on the spaces
where the function q is approximated by
For other details we refer to [14] or [11].
Further denote by the number of elements, by the number
of interior inter-element faces and the number of faces with the prescribed Neumann
boundary conditions in the discretization by
j@
j. Let e
be some numbered ordering of the set of prismatic elements and f k ,
NNC, be the ordering of the set of non-Dirichlet faces from h . For every element
we denote by NIF e i the number of interior inter-element faces and by NNC e i the
number of faces with Neumann boundary conditions imposed on the element e i . Let the
nite-dimensional space RT 0
spanned by linearly independent
basis functions v NA from the denition (1.6); let the space M 0
spanned by NE linearly independent basis functions  nally the
space M 0
spanned by NIF linearly independent basis functions  k ,
+NNC. From this Raviart-Thomas approximation we obtain the system
of linear algebraic equations in the formB @
C A =B @
where are unknowns, the
symmetric positive denite matrix block A 2 R NA;NA is given by the terms (Av
the outdiagonal block B 2 R NA;NE by (r  and the block C 2 R NA;NIF+NNC
by <
. Here n k is the outward normal vector to with respect to the face
(see [11] and [14]).
Let us denote the system matrix in (1.12) by A. The symmetric matrix A is indenite
due to the zero diagonal block of dimension NNC. The structure
of nonzero elements in the matrix from a small model problem can be seen in Figure
1. Partition the submatrix C in A as (C 1 corresponds to the
interior inter-element faces in the discretized domain and C 2 2 R NA;NNC is the face-
condition incidence matrix corresponding to the element faces with Neumann boundary
conditions. Note that every column of C 1 contains only two nonzero entries equal to 1.
The singular values of C 1 are all equal to
2 and the matrix block C 2 has orthogonal
columns. Moreover, the whole matrix block C has also singular values equal to
2 or
1. The matrix B has a special structure. The nonzero elements correspond to the face-
element incidence matrix with values equal to -1. Thus all singular values of the matrix
are equal to
5 (the matrix B is, up to the normalization coe-cients, orthogonal).
It is easy to see from the denition of approximation spaces (see [14] or [15]) that
the symmetric positive denite block is 5  5 block diagonal and it was shown in [15]
Fig. 1. Structural pattern of the matrix obtained from mixed hybrid nite element approximation of
a model problem with (to be discussed in Section 5).
that the spectrum of the matrix block A satises
are positive constants independent of the discretization parameters and
dependent on the domain and the tensor A. It is also easy to see that the system matrix
A in (1.12) is non-singular if and only if the block (B C) has a full column rank. Clearly,
if the condition
@
holds (all boundary conditions are Neumann conditions), then
the matrix block (B C) is singular, due to the fact that all sums of row elements are zero.
In other words, the function p is unique up to a constant function in the case
@
Assuming
@
it follows from the analysis presented in [15] that there exist positive
constants c 3 and c 4 such that for the singular values of the matrix block (B C) we have
Moreover, for eigenvalues of the whole symmetric indenite matrix A it follows asymptotically
are positive constants independent of system parameters.
In this paper, for solving the symmetric indenite systems (1.12), the successive reduction
to Schur complement systems is proposed. We consider three successive Schur
complement systems arising during the block elimination of unknowns which correspond
to matrix blocks A, B and C 2 respectively, or in other words, which correspond to the
elimination of the velocity variables u, the pressure variables p and of a part of the
Lagrange multipliers . While the concept of reduction to the rst and second Schur
complement systems is well known as a static condensation (described e.g. in [4], Section
V. or in [11])), the proposed reduction to the third Schur complement system seems to
be new. The main contribution of the paper consists in a detailed investigation of the
structure of nonzero entries and the spectral properties of the Schur complement ma-
trices. This enables thorough complexity analysis of the direct or iterative solution of
corresponding Schur complement systems. A brief analysis of the structure of the rst
Schur complement matrix can be found in [4] as well as a straightforward observation
that its principal leading block is a diagonal. Here we extend this analysis and discuss
the mutual relation between the number of nonzero entries in the rst Schur complement
matrix and the number of nonzeros in the system matrix (1.12). We show further that
no ll-in occurs during the process of reduction to the second and third Schur complement
system. Moreover, we prove that the number of nonzeros in both these two Schur
complements is always less than the number of nonzeros in (1.12). It is shown also that
the spectral properties of matrices in such Schur complement systems do not deteriorate
during the successive elimination. Thus an approach based on the block reduction and
subsequent iterative solution is well justied.
The outline of the paper is as follows. In Section 2, we examine the structural pattern
of resulting Schur complement matrices and give estimates for their number of nonzero
elements in terms of the discretization parameters listed above. Section 3 is devoted to the
solution of the whole indenite system (1.12) via three Schur complement reductions and
subsequent direct solution. Using the graph theoretical results we give the asymptotic
bound of the computational complexity for the Cholesky decomposition method applied
to the third Schur complement system. In Section 4, we concentrate on the spectral
properties of the Schur complement system matrices. The theoretical convergence rate of
the iterative conjugate gradient-type method in terms of the discretization parameters is
estimated. The asymptotic bounds for the computational work of the iterative solution
are given. Section 5 contains some numerical experiments illustrating the previously
developed theoretical results. Finally, we give some concluding remarks and mention
some open questions for our future work.
2. Structural properties of the Schur complement matrices. In this section
we take a closer look to the discretized indenite system and corresponding Schur complements
and we extend the brief analysis from [4]. There are several possibilities for
the choice of a block ordering in the consecutive elimination. We shall concentrate on
the block ordering which seems to be the most natural and e-cient from the point of
view of solving the nal Schur complement system by a direct solver or by a conjugate
gradient-type method. The same ordering for the elimination of the rst two blocks was
used also in [4], p. 178-181 or in [11]. Note that the static condensation is not the only
way to form the successive Schur complements. E.g., in [17] the case of the Raviart-
Thomas discretization for the closely related nodal methods was studied and reduction
to a dierent second Schur complement system was discussed.
The following simple result gives the number of nonzero elements in the triangular
part of the matrix A. By the triangular part of a matrix M we mean its upper (lower)
strict triangle diagonal. We will deal only with the structural nonzero elements
we do not take into account accidental cancellations and possible initial zero values of the
Fig. 2. Structural pattern of the Schur complement matrix A=A for A from Figure 1.
tensor of hydraulic permeability. By the structure of a matrix M we mean
0g.
Lemma 2.1. The number of nonzeros in the triangular part of A is given by
Proof. The triangular part of A has 15NE nonzeros, the block B contributes by 5NE
nonzeros, C 1 has 2NIF nonzeros and C 2 contains NNC nonzeros.
The symmetric positive denite matrix block A in (1.12) is block-diagonal, each 5  5
block corresponds to certain element in the discretization of the domain. Therefore it
is straightforward to eliminate the velocity variables u and to obtain the rst Schur
complement system with the matrix
A=A =B @
A 11 A 12 A 13
A T
A T
The structure of the matrix A=A for our example problem is shown in Figure 2. For
details we also refer to [4], p. 180-181 or [11]. For the number of nonzeros in the matrix
A=A we can show the following result.
Lemma 2.2. The number of nonzeros in the triangular part of the Schur complement
matrix A=A is equal to
Proof. Clearly, . Note that the
ll-in for A 1 C 1 is considerably higher (it is equal to 10NIF ). Further, jA 13
and
The number of
nonzeros in A 23 is equal to
Finally, note that
Observe that the directed graph of the matrix B T C 1 has the set of arcs
is an interior face of ig:
The undirected graph of C T
adjacency relation
based on the connectivity through the interior faces inside the domain. It follows
that
where e(f) and
e(f) are the two elements from E h such that
considering the relation
NIF: Putting all the partial sums together
we get the desired result.
Consider now the second Schur complement matrix

A

A 22 A 23
A T

The structure of ( A=A)=A 11 for our example matrix is shown in Figure 3. The matrix
block A 11 in the rst Schur complement matrix A=A is diagonal [4], [11]. The following
result shows that it is worth to form the Schur complement matrix ( A=A)=A 11 from
the matrix A=A since no further ll-in appears during the elimination of the block A 11
corresponding to pressure variables p and so we can further reduce the dimension of the
system.
Theorem 2.1.


A 22 A 23
A T
Proof. We have the following structural equivalences:
Fig. 3. Structural pattern of the Schur complement matrix ( A=A)=A11 for A from Figure 1.
From previous Theorem it is also easy to see that right lower block B 22 is block-diagonal
with blocks of varying size (depending on number of faces with Neumann conditions in
each element) each corresponding to a certain element in the discretization. So in the
following we will consider the third Schur complement matrix
induced by the block B 22 in the matrix ( A=A)=A 11 : We can prove a similar result
to the one given in Theorem 2.1. Therefore, the Schur complement system with the
matrix (A=A)=A 11 can be reduced to the Schur complement matrix (( A=A)=A 11 )=B 22
of dimension equal to NIF , without inducing any additional ll-in. Moreover, this can
be done using incomplete factorization procedures.
Theorem 2.2.
Proof. Using Theorem 2.1 we get
only in the
trivial singular case with jE h we get the desired result Struct(B 11
The following simple corollary gives the number of nonzero elements in the second and
third Schur complement matrices ( A=A)=A 11 and (( A=A)=A 11 )=B 22 . We shall use
these results later.
Corollary 2.1. The number of nonzeros in the triangular part of ( A=A)=A 11 is
given by
and the number of nonzeros in the triangular part of ( A=A)=A 11 )=B 22 is given by
Apart from explicit assembly of the Schur complement matrices or using them implicitly
there is another possibility which may be considered { keeping the Schur complements
in factorized form. Consider the following decomposition:
In contrast to the previous case, where the local numbering of the
faces corresponding to the individual elements did not play a role, this is not the case
now.
Theorem 2.3. Assume that all the elements within the diagonal blocks of the matrix
A are nonzero. The ll-in in
if the faces with Dirichlet
boundary conditions are numbered rst in the local ordering of each nite element.
Proof. Because of the block structure of A we can consider the individual nite
elements independently. The minimum value of the nonzero count of ^
subsequent
rows which correspond to the same nite element is 1P
it is easily checked to be minimal in this case.
Therefore, from now we assume that within each element we have rst numbered the faces
corresponding to Dirichlet boundary conditions, then the interior inter-element faces and
nally the faces with Neumann boundary conditions. The matrix (2.19) can be written
in the form
It is clear that it is more advantageous to keep most of the blocks of (2.20) in the
explicit form multiplying the factors directly. A typical example is the block ^
B,
which is a diagonal matrix. The main question here is whether we can reduce the system
further as in the previous case and at the same time keep the matrix blocks in a
factorized form. Unfortunately, there is one basic obstacle. Whereas we are able to
embed the structure of A T
into the structure of A 22 we cannot in general express
in the factorized form as
is factor which can be easily computed.
We have considered the partially factorized structure (2.20) since it is important from
a computational point of view. Using a structural prediction based on such factors is exactly
the way how to obtain the sparsity structure of explicit Schur complement matrices
22 in an e-cient way. In our implementations
we used a procedure similar to the one from [16] to get these structures.
3. Direct solution of the Schur complement systems. In the following we will
discuss the direct solution of the Schur complement systems. Namely, we will concentrate
on the system with the matrix (( A=A)=A 11 )=B 22 2 R NIF;NIF : The following theorem
gives a bound on the asymptotic work necessary to solve the linear system (1.12), which
is dominated by the decomposition of the matrix (( A=A)=A 11 )=B 22 :
Theorem 3.1. The number of arithmetic operations to solve the symmetric inde-
nite system (1.12) directly via three consecutive block eliminations and using the Cholesky
decomposition is O(NIF 2
Proof. We will only give a sketch of the proof here. The work is dominated by the
decomposition of B
12 , which has the same nonzero structure as A 22 :
Our uniform regular nite element mesh is a well-shaped mesh in a suitable sense
(see [19]). The proof of Lemma 2.2 and the statements of Theorem 2.1 and 2.2 imply
that the graph G of (( A=A)=A 11 )=B 22 is also the graph of a well-shaped mesh. Namely,
it is a bounded-degree subgraph of some overlap graph (see [18], [19]). It was shown in
[25] that the upper bound on the second-smallest eigenvalue of the Laplacian matrix of
G (the Fiedler value) is of the order O(1=NIF 2=3 using the techniques from [25]
we obtain that there exists a O(NIF 2=3 )-size bisector of G.
Therefore, G satises the so-called NIF 2=3 -separator theorem: there exist constants
such that the vertices of G can be partitioned into sets GA ; GB and
the vertex separator GC such that jGA j; jGB j  NIF and jGC j  NIF 2=3 : Moreover,
any subgraph of G satises the NIF 2=3 -separator theorem. The technique of recursive
partitioning of G called generalized nested dissection and used to reorder the considered
Schur complement matrix provides an elimination ordering with an O(NIF 2 )-bound on
the arithmetic work of Cholesky decomposition (see Theorem 6 in [12]).
Note that the explicit computation of the matrix (( A=A)=A 11 )=B 22 is necessary in
the framework of direct methods. Theorem 3.1 provides a theoretical result which is based
on spectral partitioning methods. The reordering algorithms based on the separators
obtained by the spectral partitioning techniques and applied recursively within the nested
dissection need not necessarily be the best practical approach to get a reasonable matrix
reordering. Nevertheless, experimental results with various partitioning schemes show
that high quality reorderings can be e-ciently computed in this way (see [7]). Also some
other reorderings which combine global procedures (partitioning of large meshes) and
local algorithms (like MMD) can provide reasonable strategies to nd a ll-in minimizing
permutation.
4. The conjugate gradient method applied to the Schur complement sys-
tems. In this section we concentrate on the iterative solution of the Schur complement
systems discussed in Section 3. We consider the conjugate gradient method applied
to the symmetric positive denite systems with matrices A=A, (( A=A)=A 11 ) and
22 . It is well known that the convergence rate of the conjugate gradient
method can be bounded in terms of the condition number of the corresponding Schur
complement matrix [9], [6], [26]. We show that the condition number of the matrix A=A
is asymptotically the same as the conditioning of the negative part of spectrum of the
whole indenite matrix A. Moreover, we prove that condition numbers of the matrices
grow like 1=h 2 with respect to the discretization
parameter h and they do not deteriorate during the successive eliminations. Based on
these results we estimate the number of iteration steps necessary to achieve the prescribed
tolerance in error norm reduction. We show that the number of iteration steps
necessary to reduce the error norm by the factor of " grows asymptotically like 1=h for
all three Schur complement systems. Therefore, the total number of
ops in the iterative
algorithm can be signicantly reduced due to decrease of the matrix order during the
elimination. First, we consider the following theorem.
Theorem 4.1. Let  1   be the eigenvalues of the positive
denite block A 2 R NA;NA ,  1   be the singular values of the
matrix block (B C) 2 R NA;NBC . Then for the eigenvalues of the Schur complement
Moreover, for the eigenvalues of the positive denite matrix blocks

A 22 A 23
A T

A 22 A 23
A T
The condition number of the Schur complement system matrix A=A then can be bounded
by the expression
(4.
Proof. The positive denite matrix A 1 has the spectrum 0 < 1= 1  1=
1=NA . The rst inclusion in the theorem follows from the following two inequalities 1
NA ((B C)x; (B C)x);
Similarly, from the inequalities 1
we obtain the second inclusion. The third part of the proof is completely analogous to
the second part.
Corollary 4.1. There exist positive constants c 9 and c 10 such that for the spectrum
of the Schur complement matrix A=A we have
The condition number of the matrix A=A can be
bounded as
The Schur complement system with positive denite matrix A=A can be solved iteratively
by the conjugate gradient method [9] or the conjugate residual method [6]. It
is well known that the conjugate gradient method generates the approximate solutions
which minimize the energy norm of the error at each iteration step [26], [6]. The closely
related conjugate residual method that dier only in the denition of innerproduct, on
the other hand, generates the approximate solutions which minimize their residual norm
at every iteration [6]. It is also well known fact that there exists so-called peak/plateau
connection between these methods [5] showing that there is no signicant dierence in
the convergence rate of these methods when measured by the residual norm of an approximate
solution. In our paper we use the conjugate gradient method together with
the minimal residual smoothing procedure applied on its top to get monotonic residual
norms [28]. Applying such technique allows better monitoring of the convergence
by residual norm and it is mathematically equivalent to the residual minimizing conjugate
residual method [6]. The computational cost of this technique is minimal and it
costs only two inner products and one vector update per iteration. In the framework
of iterative methods the number of operations in matrix-vector products is what is usually
the most important. These products, performed repeatedly in each iteration loop,
contribute in a substantial way to the nal e-ciency of iterative solver. When solving
the system with Schur complement matrix A=A the number of
ops per iteration for
an unpreconditioned method is dominated by the matrix vector multiplication with the
matrix A=A. Its number of nonzeros was given by Lemma 2.2. Moreover, using the
estimates (1.13) and (1.14), the condition number of the Schur complement matrix A=A
can be bounded by the term O( 3
Consequently, the number of
ops for conjugate
gradients necessary to achieve a reduction by " is of order
Assuming overestimates we obtain the asymptotic
estimate of order O(NE 3
The previous considerations did not take into account the Schur complement systems
with matrices (( A=A)=A 11 ) and ((( A=A)=A 11 )=B 22 ). The convergence rate of the
iterative conjugate gradient method applied to the second and third Schur complement
systems depend analogously on the condition number of the Schur complement matrices
[9], [6], [26]. The analysis of the spectrum of the matrix (( A=A)=A 11 ) is given in the
following theorem.
Theorem 4.2. Let  1   be the eigenvalues of the positive
denite block A 2 R NA;NA ,  1   be the singular values of the
matrix block (B C) 2 R NA;NBC . Then for the spectrum of the Schur complement matrix
A=A)=A 11 we have
Consequently, the condition number of the matrix ( A=A)=A 11 can be bounded as follows
Proof. From the denition of the Schur complement matrix ( A=A)=A 11 and the
statement of Theorem 4.1 we have

A 22 A 23
A T
The bound for the minimal eigenvalue can be obtained considering the following result
(see [20], p.201):
A 11 A 12 A 13
A T
A T

A

A
Then from the interlacing property of the eigenvalue set of symmetric matrix A=A (see
e.g. [8]) it follows




Considering the previous inequalities we get the lower bound for the minimal eigenvalue
of the matrix ( A=A)=A 11 , which completes the proof.
We have shown that the condition number of the Schur complement system matrix
A=A)=A 11 is bounded by a multiple of the condition number of the matrix A=A.
Therefore the number of iteration steps for the conjugate gradient method necessary to
reduce the error norm(or after smoothing the residual norm) by some factor is asymptotically
the same as before. The complexity of the matrix-vector multiplication is lower
and according to Corollary 2.1 is of the order
Assuming again the overestimates (NIF we obtain the
asymptotic estimate O(NE). The total number of
ops for the conjugate gradients
or the conjugate residual method necessary to achieve a reduction by the factor " is
then again of order O(NE 3
NE). From the statements of Theorem 4.1 and Theorem
4.2 it is clear that the reduction to the Schur complement systems does not aect the
asymptotic conditioning of the positive denite matrices A=A and ( A=A)=A 11 : The
same is true for the spectral properties of the third Schur complement system with
the matrix (( A=A)=A 11 )=B 22 : Since the proof is completely analogous to the proof of
Theorem 4.2 we shall present only the following statement (cf. [10], p. 256).
Theorem 4.3. The condition number of (( A=A)=A 11 )=B 22 is bounded by the
condition number of the matrix ( A=A)=A 11
In the following we present two additional results concerning the the matrix-vector
multiplications with Schur complement matrices. Theorem 4.4 compares the number of
nonzeros in the Schur complement matrices ( A=A)=A 11 and (( A=A)=A 11 )=B 22 to the
number of nonzeros in the original matrix A:
Theorem 4.4. The number of nonzero entries in the matrix (( A=A)=A 11 ) or the
matrix ((( A=A)=A 11 )=B 22 ) is smaller than the number of nonzeros in the matrix A:
Proof. Using the fact that
it follows from Lemma 2.1 and Theorem 2.1 that
Clearly, the number of nonzeros in the matrix ((( A=A)=A 11 )=B 22 ) is even smaller.
Note that the number of nonzeros in the original matrix A can be smaller or larger
than the corresponding number of nonzeros in the matrix A=A: Consider now the factorized
Schur complement in the form (2.20). It can be shown also that there is no clear
winner between the number of
oating-point operations to multiply a dense vector by the
or the number of operations to get a product of
a matrix (( A=A)=A 11 )=B 22 with a dense vector of appropriate dimension, respectively.
The result depends on the shape of the domain and its boundary conditions. Never-
theless, the following Theorem 4.5 shows that if we do not form the Schur complement
explicitly it is worth to use the factorized form (2.19) and the reordering of the Schur
complement from Theorem 2.3 instead of its implicit form.
Theorem 4.5. Let v be a dense vector. The number of
oating-point operations
to compute
smaller than the number of
oating-point
operations to computeB @
Proof. Taking into account the local ordering from Theorem 2.3 the dierence between
these two quantities can be bounded by
2NIF NNC  0:
5. Numerical experiments. In the following we present numerical experiments
which illustrate the results developed in the theoretical part of the paper.
Two model potential
ow problems (1.1) and (1.2) in a rectangular domain with
Neumann conditions prescribed on the bottom and on the top of the domain have been
considered. Dirichlet conditions that preserve the nonsingularity of the whole system
matrix A were imposed on the rest of the boundary. The choice of boundary conditions
in these examples is motivated by our application and it comes from a modelling of a
conned aquifer (see [3]) between two impermeable layers.
In order to verify the theoretical results derived in previous sections we will restrict
our attention rst to the simplest geometrical shape - cubic domain and report the results
obtained from a uniformly regular mesh renement. In practical situations, however, relatively
thin aquifers with possible cracks in the rock are frequently modelled, and so the
number of Neumann conditions may represent a big portion of the whole boundary. As
our second model example, we consider a rectangular domain discretized by 6 layers of
Model potential
uid
ow problem - cubic domain
Discretization parameters Matrix dimensions
h, NE NIF NNC A A=A ( A=A)=A 11 (( A=A)=A 11 )=B 22
1/5, 250 525 100 2125 875 625 525
1/10, 2000 4600 400 17000 7000 5000 4600
1/15, 6750 15975 900 57375 23625 16875 15975
1/20, 16000 38400 1600 136000 56000 40000 38400
1/30, 54000 131400 3600 459000 189000 135000 131400
1/35, 87750 209475 4900 728875 300125 214375 209475
1/40, 128000 313600 6400 1088000 448000 320000 313600

Table
Model potential
uid
ow problem - realistic domain
Discretization parameter Matrix dimension
95x95x6 251560 36100 937460 395960 287660 251560
elements in the mesh. As we will see later, the reduction to the third Schur complement
proposed in this paper can become even more signicant than for the cubic domain.
Prismatic discretizations of domains with NE elements were used [14], [11]. For the
cubic domain we have then Discretization parameters h, NE, NIF , NNC,
dimension N of the resulting indenite system matrix A and the dimensions of the corresponding
Schur complement matrices A=A, ( A=A)=A 11 and (( A=A)=A 11 )=B 22 are
given in Table 1 for a cubic domain and in Table 2 for a more realistic domain. We note
again that the dierence between dimensions of the second and third Schur complement
matrix is signicantly larger in the case of modelling of thin layers that arise regularly
in our application.
For the example of a cubic domain the spectral properties of the matrix blocks A
and (B C) as well as of the whole symmetric indenite matrix A have been investigated.
The extremal positive and negative eigenvalues of the matrix A and the extremal singular
values of the block (B C) (squared roots of the extremal eigenvalues of the matrix
were approximated by a reduction to the symmetric tridiagonal form of
the matrix using 1500 steps of the symmetric Lanczos algorithm [8] and by a subsequent
Spectral properties of the system matrix and its blocks - problem with a cubic domain
matrix blocks spectral properties eigenvalues of the matrix A
NE spectrum of A sing. values of (BjC) negative part positive part
2000 [0.33e-2, 0.2e-1] [0.927e-1, 2.64] [-2.64, -0.898e-1] [0.335e-2, 2.64]
16000 [0.66e-2, 0.4e-1] [0.467e-1, 2.64] [-2.64, -0.413e-1] [0.679e-2, 2.65]
54000 [0.99e-2, 0.6e-1] [0.312e-1, 2.65] [-2.64, -0.241e-1] [0.104e-1, 2.65]
128000 [0.13e-1, 0.8e-1] [0.234e-1, 2.65] [-2.64, -0.152e-1] [0.136e-1, 2.65]
eigenvalue computation of the resulting tridiagonal matrix using the LAPACK double
precision subroutine DSYEV [1]. Extremal eigenvalues of the diagonal matrix block A
were computed directly by the LAPACK eigenvalue solver element by element. It can
be seen that the computed extremal eigenvalues of the block A are in perfect agreement
with the theory (see Table 3). Similarly, we can observe approximately a linear decrease
of the computed minimal singular value of the matrix block (B C) with respect to the
mesh discretization parameter h. From the computed extremal eigenvalues of the whole
indenite system A we can conclude that even if our mesh size parameters h are rather
small and give rise to very large system dimensions (see Table 1), they are outside of
the asymptotic inclusion set (1.15). Indeed for our example and our mesh size interval
we have c 1 =h  c 4 , c 2 =h  c 4 and with the exception of
using Lemma 2.1 in [22], pp. 3-4 (see also [15]) we obtain the inclusion
set in the form
which is in good agreement with the results in Table 3.
Using the same technique we have approximated the extremal eigenvalues of the
Schur complement matrices A=A, ( A=A)=A 11 and (( A=A)=A 11 )=B 22 coming from
a problem on a cubic domain. From Table 4 it can be seen that the inclusion set for
the extremal eigenvalues of the rst Schur complement matrix A=A coincides with the
bounds given in Theorem 4.1. We can see that the extremal eigenvalues of the second
Schur complement matrix ( A=A)=A 11 are bounded by the extremal eigenvalues of the
matrix A=A. Similarly, the extremal eigenvalues of the third Schur complement matrix
are bounded by the extremal eigenvalues of the matrix ( A=A)=A 11 .
This behaviour is in accordance with the asymptotic bounds given in Theorem 4.2 and
Theorem 4.3.
The smoothed conjugate gradient method has been applied to the resulting three
Schur complement systems (see also the discussion in previous section). Unpreconditioned
and also preconditioned versions with the IC(0) preconditioner [23], [24] for the
solution of these symmetric positive denite systems have been used. For the solution of
Spectral properties of Schur complement matrices - problem with a cubic domain
spectral properties of Schur complement matrices
2000 [0.182e1, 0.173e4] [0.251e1, 0.596e3] [0.272e1, 0.596e3]
54000 [0.693e-1, 0.579e3] [0.966e-1, 0.199e3] [0.992e-1, 0.199e3]
128000 [0.293e-1, 0.434e3] [0.409e-1, 0.149e3] [0.417e-1, 0.149e3]
the whole indenite system the minimal residual method has been used. For the preconditioned
version the positive denite block-diagonal preconditioning with ILUT(0,20) for
the decomposition of the block corresponding to constraints (see e.g. [22], [21]) was used.
The choice of ILUT(0,20) was motivated by our eort to obtain rather precise factorization
with restricted memory requirements which should be close to the full decomposition
of the block (B C) T (B C). This preconditioner was found generally better than the in-
denite block-diagonal preconditioning with the same ILUT(0,20) decomposition or than
the indenite preconditioner discussed in [13] or [21]. The initial approximation x 0 was
set to zero, the relative residual norm krnk
used as the stopping criterion.
For the implementation details of iterative solvers we refer to [6]. Our experiments were
performed on an SGI Origin 200 with processor R10000. In Table 5 and Table 6 we
consider iteration counts and CPU times in the minimal residual method (unprecondi-
tioned/preconditioned) applied to the whole system (1.12) and in the conjugate gradient
method (unpreconditioned/preconditioned) applied to the Schur complement systems
with the matrices A=A, ( A=A)=A 11 and (( A=A)=A 11 )=B 22 for a model problem with
a cubic and more realistic domain, respectively. The dependence of the iteration counts
presented in all columns of Table 5 corresponds surprisingly well to the theoretical order
O( 3
NE). The convergence behaviour of the smoothed conjugate gradient method applied
to the third Schur complement system with the matrix (( A=A)=A 11 )=B 22 for this
case is presented in Figure 4. From the results in Table 5 and Table 6 it follows that while
the gain from the solution of the third Schur complement system is rather moderate in
the case of a cubic domain and in the case of the realistic
at domain it becomes more
signicant.
6. Conclusions. Successive block Schur complement reduction for the solution of
symmetric indenite systems has been considered in the paper. It was shown that due
to the particular structure of matrices which arise from mixed-hybrid nite element discretization
of the potential
uid
ow problem, the resulting Schur complement matrices
remain sparse. Moreover, their spectral properties do not deteriorate and the iterative
conjugate gradient method can be successfully applied. Theoretical bounds for the
Number of iterations of the conjugate gradient method - problem with a cubic domain
unpreconditioned/preconditioned CG applied to matrix
NE A A=A ( A=A)=A 11 (( A=A)=A 11 )=B 22
2000 608/76 154/35 87/32 80/32
6.43/1.56 0.76/0.25 0.30/0.16 0.25/0.14
48.17/11.91 3.86/1.50 1.51/0.92 1.30/1.01
16000 1031/138 288/67 164/63 155/63
54000 1358/188 418/95 234/93 228/93
926.98/218.78 104.76/36.92 39.88/24.04 37.94/23.13
128000 1637/229 546/122 303/122 298/122

Table
Number of iterations of the conjugate gradient method - realistic model example
unpreconditioned/preconditioned CG applied to matrix
NE A A=A ( A=A)=A 11 (( A=A)=A 11 )=B 22
50700 2053/336 810/172 448/166 421/166
86700 2959/403 1042/222 578/214 543/214
132300 3420/447 1272/271 706/262 663/262
iteration number
relative
residual
norms
unpreconditioned and smoothed conjugate gradient method applied to ((-A/A)/A 11 )/B 22
Fig. 4. Convergence of the smoothed conjugate gradient method applied to the third Schur complement
system
convergence rate of this method in terms of the discretization parameters have been developed
and tested on a model problem example. Numerical experiments indicate that
the given theoretical bounds on the eigenvalue set are realistic not only for the system
matrix and its blocks, but also for the Schur complement matrices. The iteration counts
for the conjugate gradient method are also in a good agreement with the theoretical
predictions. Direct solution of the third Schur complement system is also a possible al-
ternative. Nevertheless, its comparison with iterative solvers is outside the scope of this
paper.
In case of structured grids, a geometric multigrid solver and/or preconditioner for
solving the nal Schur complement system can be used. Namely, the stencil from the rst
Schur complement which expresses element-element connectivity in the domain (see proof
of Lemma 2.2) remains unchanged after the subsequent two reduction and an appropriate
method could be based on that.
Another approach for the solution of symmetric indenite systems seems to be
promising. As was pointed out in [2], the classical null-space algorithm can be imple-
mented. QR factorization of the o-diagonal block (B C) is considered and the solution
of the indenite system is transformed to the solution of a block lower triangular sys-
tem, where the subproblem corresponding to the diagonal block can be solved using the
factorization or an iterative conjugate gradient-type algorithm. This approach
has the advantage of performing the matrix-vector multiplication by the Q factor using
elementary Householder transformations. Although the Q factor may be structurally
full, the elementary Householder vectors may be quite sparse. Moreover, a roundo error
analysis of the algorithm can be carried out.
7.

Acknowledgment

. Authors would like to thank Michele Benzi for careful reading
of manuscript and anonymous referees for their many useful comments which significantly
improved the presentation of the paper. We are indebted to Jir Muzak from the
Department of Mathematical Modelling in DIAMO, s.e., Straz pod Ralskem for providing
us with a model numerical example for the experimental part of this paper and to Jorg
Liesen for giving us the reference [20]. This work was supported by the Grant Agency of
the Czech Republic under grant 201/98/P108 and by the grant AS CR A2030706.



--R

LAPACK User's Guide SIAM
The use of QR factorization in sparse quadratic programming.
Dynamics of Fluids in Porous Media.
Mixed and Hybrid Finite Element Methods.
Relations between Galerkin and norm-minimizing iterative methods for solving linear systems
Iterative methods for large
Geometric mesh partitioning: implementation and exper- iments
Matrix Computations.
Method of conjugate gradients for solving linear systems.
Accuracy and Stability of Numerical Algorithms.

Generalized nested dissection



Sparse QR factorization with applications to linear least squares problems.
Approximate Schur complement preconditioning of the lowest-order nodal discretizations
Automatic mesh partitioning
Geometric separators for
Schur Complement and Statistics.

A preconditioned iterative method for saddle point problems.
Iterative Methods for Sparse Linear Systems.
ILUT: A dual threshold incomplete ILU factorization.
Spectral partitioning works: Planar graphs and
Parallel iterative solution methods for linear systems arising from discretized PDE's.


--TR
indefinite linear systems;preconditioned conjugate residuals;potential fluid flow problem;sparse linear systems;finite element matrices