Graphics Reference
In-Depth Information
The main reason that we don't use this method is that although
A
has
very few non-zeros,
L
can have a lot. In three dimensions the amount of
fill-in
(extra non-zeros in
L
) is particularly bad;
direct solvers
that take
this approach generally fail on 3D problems from lack of memory.
Incomplete Cholesky tackles this problem with a very simple idea:
whenever the Cholesky algorithm tries to create a new non-zero in a lo-
cation that is zero in
A
, cancel it—keep it as zero. On the one hand, the
resulting
L
is just as sparse as
A
, and memory is no longer an issue. On
the other hand, we deliberately made errors:
A
=
LL
T
now.
However
,
hopefully the incomplete factorization is close enough to
A
that doing the
solves in Equation (4.13) is close enough to applying
A
−
1
so that we have
a useful preconditioner for PCG!
Technically, performing incomplete Cholesky only allowing non-zeros
in
L
where there are non-zeros in
A
is called level zero: IC(0). There are
variations that allow a limited number of non-zeros in other locations, but
we will not broach that topic here.
To make this more precise, IC(0) constructs a lower-triangular matrix
L
with the same non-zero pattern as the lower triangle of
A
, such that
LL
T
=
A
wherever
A
is non-zero. The only error is that
LL
T
is non-zero
in some other locations where
A
is zero.
Assume we order our grid cells (and the corresponding rows and
columns of
A
) lexicographically, say along the
i
-dimension first, then the
j
-dimension, and finally the
k
-dimension.
5
Suppose we split
A
up into its
strict lower triangle
F
and diagonal
D
:
A
=
F
+
D
+
F
T
.
Then, it can be shown for the particular
A
we're solving—though we won't
show it here—that the IC(0) factor
L
is of the form
L
=
FE
−
1
+
E,
where
E
is a diagonal matrix. That is, all we need to compute and store
are the diagonal entries of
L
, and we can infer the others just from
A
!
Crunching through the algebra gives the following formulas for comput-
ing the entries along the diagonal of
E
. In two dimensions,
E
(
i,j
)
=
A
(
i,j
)
,
(
i,j
)
−
(
A
(
i−
1
,j
)
,
(
i,j
)
/E
(
i−
1
,j
)
)
2
−
(
A
(
i,j−
1)
,
(
i,j
)
/E
(
i,j−
1)
)
2
.
5
It turns out that the order in which we take the dimensions doesn't actually change
anything for IC(0) applied to our particular matrix.
