Graphics Reference
In-Depth Information
floating-point numbers, 64-bit double precision is
strongly
recommended for
PCG—at the very least, for the scalars (in particular when accumulating
the dot-products
4
). Single precision rounding errors can cause a significant
slow-down in convergence.
4.3.3
Incomplete Cholesky
We still have the question of defining the preconditioner. From the stand-
point of convergence the perfect preconditioner would be
A
−
1
, except that's
obviously far too expensive to compute. The true ideal is something that is
both fast to compute and apply and is effective in speeding up convergence,
so as to minimize the total solution time.
There are many, many choices of preconditioner, with more being in-
vented each year. Our default choice, though, is quite an old preconditioner
from the
incomplete Cholesky
(IC) family. It's both simple to implement
and fairly ecient, and it is quite robust in handling irregular domains
(like the shape of a liquid splash). Its chief problems are that it's hard to
parallelize effectively and that it's not optimally scalable (the number of
iterations required for PCG steadily increases with grid size); more com-
plex approaches such as
domain decomposition
or
multigrid
may become
useful if you want to scale up your simulator to run on large problems on
multi-core machines or clusters, but these methods lie beyond the scope of
this topic.
Recall how you might directly solve a system of linear equations with
Gaussian elimination. That is, you perform row reductions on the system
until the matrix is upper-triangular, and then use back substitution to get
the solution one entry at a time. Mathematically, it turns out that this
is equivalent to factoring the matrix
A
as the product of a lower- and an
upper-triangular matrix and then solving the triangular systems one after
the other. In the case of a symmetric positive definite
A
, we can actually
do it so that the two triangular matrices are transposes of each other:
A
=
LL
T
.
This is called the
Cholesky
factorization. The original system
Ap
=
b
is
the same as
L
(
L
T
p
)=
b
, which we can solve as
solve
Lq
=
b
with forward substitution
,
(4.13)
L
T
p
=
q
solve
with backward substitution
.
4
If you do use single precision floating-point values for your vectors (pressure, etc.),
you should investigate the BLAS routine
dsdot
to compute the dot-product in double
precision.
