Graphics Reference
In-Depth Information
•
Set initial guess
p
= 0 and residual vector
r
=
b
(If
r
=0then
return
p
)
•
Set auxiliary vector
z
=
applyPreconditioner
(
r
), and search
vector
s
=
z
•
σ
=
dotproduct
(
z, r
)
•
Loop until done (or maximum iterations exceeded):
•
Set auxiliary vector
z
=
applyA
(
s
)
•
α
=
σ/
dotproduct
(
z, s
)
•
Update
p
←
p
+
αs
and
r
←
r
−
αz
•
If max
|
r
|≤
tol
then return
p
•
Set auxiliary vector
z
=
applyPreconditioner
(
r
)
•
σ
new
=
dotproduct
(
z, r
)
•
β
=
σ
new
/σ
•
Set search vector
s
=
z
+
βs
•
σ
=
σ
new
•
Return
p
(and report iteration limit exceeded)
Figure 4.5.
The preconditioned conjugate gradient (PCG) algorithm for solving
Ap
=
d
.
Note that PCG needs storage for an “auxiliary” vector
z
and a “search”
vector
s
(thesamesizeas
p
,
b
,and
r
), and calls subroutine
applyA
to
multiply the coecient matrix
A
times a vector and subroutine
applyPreconditioner
to multiply
M
by a vector (which we will talk about
next). Also note that the new residual in each iteration is incrementally
updated from the previous residual—not only is this more ecient than
calculating it from scratch each time, but it in fact tends to reduce the
number of iterations required due to some interesting round-off error in-
teractions. Since PCG tends to be the most time-consuming part of the
simulation, it pays to use optimized
BLAS
3
routines for the vector opera-
tions here. Finally, it's worth noting that although most of the rest of the
fluid simulator may be effectively implemented using 32-bit single precision
3
The BLAS, or basic linear algebra subroutines, provide a standardized API for sim-
ple operations like dot-products, scaling and adding vectors, multiplying dense matrices,
etc. Every significant platform has at least one implementation, which generally is in-
tensely optimized, making full use of vector units, cache prefetching, multiple cores, and
the like. It can be very dicult (if not impossible without using assembly language)
to match the eciency attained by the BLAS for vectors of any appreciable size, so it
generally pays to exploit it.
