White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI
Подождите немного. Документ загружается.
8.6. PARALLEL ADI SCHEMES 341
where
P
1
(L + Y )
1
[L + (L K)(L + K)
1
] and
P
1
Q = (L + Y )
1
(L K)(L + K)
1
(L Y )=
Thus, the convergence can be analyzed by either requiring some norm of P
1
Q
to be less than one, or by requiring D to be SPD and D = P Q to be a P-
regular splitting. The following theorem uses the P-regular splitting approach.
Theorem 8.6.1 (ADI Splitting Convergence) Consider the ADI algorithm
where
is some positive constant. If D> K and Y are SPD and is such
that
L +
1
2
(Y K + KY ) is positive definite, then the ADI splitting is P-regular
and must converge to the solution.
Proof. Use the above splitting associated with the ADI algorithm
P
1
= (L + Y )
1
[L + (L K)(L + K)
1
]
= (
L + Y )
1
[(L + K) + (L K)](L + K)
1
= (L + Y )
1
2(L + K)
1
=
Thus, P =
1
2
(L + K)(L + Y ) and Q = D + P = Y K + P so that by
the symmetry of
K and Y
P
W
+ Q = (
1
2
(L + K)(L + Y ))
W
Y K +
1
2
(L + K)(L + Y )
=
1
2
(L + Y
W
)(L + K
W
) Y K +
1
2
(L + K)(L + Y )
=
L +
1
2
(Y K + KY )=
Example 2 in Section 8.2 implies that, for suitably large > P
W
+ Q =
L +
1
2
(Y K +KY ) will be positive definite. The parallel aspects of this method
are that L + K and L + Y are essentially block diagonal with tridiagonal
matrices, and therefore, the solve s teps have many independent substeps. The
ADI method may also be applied to three space dimension problems, and also
there are several domain decomposition variations of the tridiagonal algorithm.
ADI in Three Space Variables.
Consider the partial di
erential equation in three variables x
{{
x
||
x
}}
= i= Discretize this so that the algebraic problem D{ = g can be broken
into parts associated with three directions and
D = K + Y + Z where Z is
associated with the z direction. Three splittings are
D = (L + K) (L Y Z )
D = (L + Y ) (L K Z )
D = (L + Z ) (L K Y )=
© 2004 by Chapman & Hall/CRC
342 CHAPTER 8. CLASSICAL METHODS FOR AX = D
The ADI scheme has the form
for m = 0, maxm
solve (
L + K){
p+1@3
= g + (L Y Z ){
p
solve (L + Y ){
p+2@3
= g + (L K Z ){
p+1@3
solve (L + Z ){
p+1
= g + (L Y K){
p+2@3
test for convergence
endloop.
For this three variable case one can easily prove an analogue to Theorem 8.6.1.
Tridiagonal with Domain Decomposition and Interface Blo cks.
Consider the tridiagonal matrix with dimension
q = 3p + 2= Reorder the
unknowns so that unknowns for
l = p + 1 and l = 2p + 2 are listed last. The
reordered coe
!cient matrix for p = 3 will have the following nonzero pattern
5
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
7
{ {
{ { {
{ { {
{ { {
{ { {
{ { {
{ { {
{ { {
{ {
{ { {
{ { {
6
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
8
=
This the form associated with the Schur complement where E is the block
diagonal with three tridiagonal matrices. There will then be three independent
tridiagonal solves as substeps for solving this domain decomposition or "arrow"
matrix.
Tridiagonal with Domain Decomposition and No Interface Blocks.
Consider the tridiagonal matrix with dimension
q = 3p= Multiply this
matrix by the inverse of the block diagonal of the matrices D(1 : p> 1 : p),
D(p + 1 : 2p> p + 1 : 2p) and D(2p + 1 : 3p> 2p
+ 1 : 3p)= The new matrix
has the following nonzero pattern for p = 3
5
9
9
9
9
9
9
9
9
9
9
9
9
7
1 {
1 {
1 {
{
1 {
{
1 {
{
1 {
{
1
{ 1
{ 1
6
:
:
:
:
:
:
:
:
:
:
:
:
8
=
© 2004 by Chapman & Hall/CRC
8.6. PARALLEL ADI SCHEMES 343
Now reorder the unknowns so that unknowns for
l = p> p + 1> 2p and 2p + 1
are listed last. The new coe
!cient matrix has the form for p = 3
5
9
9
9
9
9
9
9
9
9
9
9
9
7
1 {
1 {
1 { {
1 {
1 {
1 {
{
1 {
{ 1 {
{
1
6
:
:
:
:
:
:
:
:
:
:
:
:
8
=
This matrix has a bottom d iagonal block, which can be permuted to a 4 × 4
tridiagonal matrix. The top diagonal block is the identity matrix with dimen-
sion 3
p 4 so that the new matrix is a 2 × 2 block upper triangular matrix.
More details on this approach can be found in the paper by N. Mattor, T. J.
Williams and D. W. Hewett [15].
8.6.1 Exercises
1. Generalize Theorem 8.6.1 to the case where is a diagonal matrix with
positive diagonal components.
2. Generalize Theorem 8.6.1 to the three space variable case.
3. Generalize the tridiagonal algorithm with s 1 interface nodes and q =
sp + s 1 where there are s blocks with p unknowns per block. Implement
this in an MPI code using s processors.
4. Generalize the tridiagonal algorithm with no interface nodes and
q = sp
where there are s blocks with p unknowns per block. Implement this in an
MPI code using s processors.
© 2004 by Chapman & Hall/CRC
Chapter 9
Krylov Methods for Ax = d
In this chapter
we show each iterate of the conjugate gradient method can be expressed as the
initial guess plus a linear combination of the
D
l
u({
0
) where u
0
= u({
0
) = gD{
0
is the initial residual and D
l
u({
0
) are called the Krylov vectors. Here {
p
= {
0
+
f
0
u
0
+ f
1
Du
0
+ · · · + f
p1
D
p1
u
0
and the coe!cients are the unknowns. They
can b e determined by either requiring M({
p
) =
1
2
({
p
)
W
D{
p
({
p
)
W
g> where D
is a SPD matrix, to be a minimum, or to require U ({
p
) = u({
p
)
W
u({
p
) to be a
minimum. The first case gives rise to the conjugate gradient method (CG), and
the second case generates the generalized minimum residual method (GMRES).
Both methods have many variations, and they can be accelerated by using
suitable preconditioners,
P, where one applies these methods to P
1
D{ =
P
1
g= A number of preconditioners will be described, and the parallel aspects
of these methods will be illustrated by MPI codes for preconditioned CG and
a restarted version of GMRES.
9.1 Conjugate Gradient Method
The conjugate gradient method as described in Sections 3.5 and 3.6 is an en-
hanced version of the method of steepest descent. First, the multiple residuals
are used so as to increase the dimensions of the underlying search set of vectors.
Second, conjugate directions are used so that the resulting algebraic systems for
the coe
!cients will be diagonal. As an additional benefit to using the conjugate
directions, all the coe
!cients are zero except for the l ast one. This means the
next iterate in the conjugate gradient method is the previous iterate plus a con-
stant times the last search direction. Thus, not all the search directions need
to be stored. This is in contrast to the GMRES method, which is applicable to
problems where the coe
!cient matrix is not SPD.
The implementation of the conjugate gradient method has the form given
below. The steepest descent in the direction
s
p
is given by using the parameter
> and the new conjugate direction s
p+1
is given by using the parameter = For
345
© 2004 by Chapman & Hall/CRC
The conjugate gradientmethodwas introduced in Chapter 3.
346 CHAPTER 9. KRYLOV METHODS FOR AX = D
each iteration there are two dot products, three vector updates and a matrix-
vector product. These substages can be done in parallel, and often one tries to
avoid storage of the full coe
!cient matrix.
Conjugate Gradient Method.
Let
{
0
be an initial guess
u
0
= g D{
0
s
0
= u
0
for m = 0, maxm
= (u
p
)
W
u
p
@(s
p
)
W
Ds
p
{
p+1
= {
p
+ s
p
u
p+1
= u
p
Ds
p
test for convergence
= (u
p+1
)
W
u
p+1
@(u
p
)
W
u
p
s
p+1
= u
p+1
+ s
p
endloop.
The connection with the Krylov vectors D
l
u({
0
) evolves from expanding the
conjugate gradient loop.
{
1
= {
0
+
0
s
0
= {
0
+
0
u
0
u
1
= u
0
0
Ds
0
= u
0
0
Du
0
s
1
= u
1
+
0
s
0
= u
1
+
0
u
0
{
2
= {
1
+
1
s
1
= {
1
+
1
¡
u
1
+
0
u
0
¢
=
{
0
+
0
u
0
+
1
¡
u
0
0
Du
0
+
0
u
0
¢
=
{
0
+ f
0
u
0
+ f
1
Du
0
.
.
.
{
p
= {
0
+ f
0
u
0
+ f
1
Du
0
+ · · · + f
p1
D
p1
u
0
=
An alternative definition of the conjugate gradient method is to choose the
co e
!cients of the Krylov vectors so as to minimize M({)
M({
p+1
) = min
f
M({
0
+ f
0
u
0
+ f
1
Du
0
+ · · · + f
p
D
p
u
0
)=
Another way to view this is to define the Krylov space as
N
p
{{ | { = f
0
u
0
+ f
1
Du
0
+ · · · + f
p1
D
p1
u
0
> f
l
5 R}=
The Krylov spaces have these very useful properties:
N
p
N
p+1
(9.1.1)
DN
p
N
p+1
and (9.1.2)
N
p
{{ | { = d
0
u
0
+ d
1
u
1
+ · · · + d
p1
u
p1
> d
l
5 R}= (9.1.3)
© 2004 by Chapman & Hall/CRC
9.1. CONJUGATE GRADIENT METHOD 347
So the alternate definition of the conjugate gradient method is to choose
{
p+1
5 {
0
+ N
p+1
so that
M({
p+1
) = min
|5{
0
+N
p+1
M(|)= (9.1.4)
This approach can be shown to be equivalent to the original version, see C. T.
In order to gain some insight to this let
{
p
5 {
0
+ N
p
b e such that it
minimizes M(|) where | 5 {
0
+ N
p
= Let | = {
p
+ w} where } 5 N
p
and use
the identity in (8.4.3)
M({
p
+ w}) = M({
p
) u({
p
)
W
w} +
1
2
(
w})
W
D(w})=
Then the derivative of M({
p
+ w}) with respect to the parameter w evaluated at
w = 0 must be zero so that
u({
p
)
W
} = 0 for all } 5 N
p
= (9.1.5)
Since N
p
contains the previous residuals,
u({
p
)
W
u({
o
) = 0 for all o ? p=
Next consider the dierence between two iterates given by (9.1.4), {
p+1
=
{
p
+ z where z 5 N
p+1
. Use (9.1.5) with p replaced by p + 1
u({
p+1
)
W
} = 0 for all } 5 N
p+1
(g D({
p
+ z))
W
} = 0
(
g D{
p
)
W
} z
W
D
W
} = 0
u({
p
)
W
} z
W
D} = 0=
Since N
p
N
p+1
> we may let } 5 N
p
and use (9.1.5) to get
z
W
D} = ({
p+1
{
p
)
W
D} = 0= (9.1.6)
This implies that z = {
p+1
{
p
5 N
p+1
can be expressed as a multiple of
u({
p
) + b} for some b} 5 N
p
= Additional inductive arguments s how for suitable
and
{
p+1
{
p
= (u({
p
) + s
p1
)= (9.1.7)
This is a noteworthy property because one does not need to store all the con-
jugate search directions
s
0
> · · · > s
p
=
The proof of (9.1.7) is done by mathematical induction. The term {
0+1
5
{
0
+ N
0+1
is such that
M({
0+1
) = min
|5{
0
+N
0+1
M(|)
= min
f
0
M({
0
+ f
0
u
0
)=
© 2004 by Chapman & Hall/CRC
Kelley [11, Chapter 2].
348 CHAPTER 9. KRYLOV METHODS FOR AX = D
Thus,
{
1
= {
0
+
0
u
0
for an appropriate
0
as given by the identity in (8.4.3).
Note
u
1
= g D{
1
= g D({
0
+
0
u
0
)
= (
g D{
0
)
0
Du
0
= u
0
0
Ds
0
=
The next step is a little more interesting. By definition {
1+1
5 {
0
+ N
1+1
is
such that
M({
1+1
) = min
|5{
0
+N
1+1
M(|)=
Properties (9.1.5) and (9.1.6) imply
(u
1
)
W
u
0
= 0
(
{
2
{
1
)
W
Du
0
= 0=
Since {
2
{
1
5 N
2
> {
2
{
1
= (u
1
+ s
0
) and
(
{
2
{
1
)
W
Du
0
= 0
(
u
1
+ s
0
)
W
Du
0
= 0
(
u
1
)
W
Du
0
+ (s
0
)
W
Du
0
= 0=
Since u
1
= u
0
0
Ds
0
and (u
1
)
W
u
0
= 0>
(u
1
)
W
u
0
= (u
0
0
Ds
0
)
W
u
0
= (u
0
)
W
u
0
0
(Ds
0
)
W
u
0
= 0=
So, if u
0
6= 0> (Ds
0
)
W
u
0
6= 0 and we may choose
=
(u
1
)
W
Du
0
(Ds
0
)
W
u
0
=
The inductive step of the formal proof is similar.
Theorem 9.1.1 (Alternate CG Method) Let
D be SPD and let {
p+1
be the
alternate definition of the conjugate gradient method in (9.1.4). If the residuals
are not zero vectors, then equations (9.1.5-9.1.7) are true.
The utility of the Krylov approach to both the conjugate gradient and the
generalized residual methods is a very nice analysis of convergence properties.
These are based on the following algebraic identities. Let
{ be the solution of
D{ = g so that u({) = g D{ = 0> and use the identity in line (8.4.3) for
symmetric matrices
M({
p+1
) M({) =
1
2
(
{
p+1
{)
W
D({
p+1
{)
=
1
2
°
°
{
p+1
{
°
°
2
D
=
© 2004 by Chapman & Hall/CRC
9.1. CONJUGATE GRADIENT METHOD 349
Now write the next iterate in terms of the Krylov vectors
{ {
p+1
= { ({
0
+ f
0
u
0
+ f
1
Du
0
+ · · · + f
p
D
p
u
0
)
=
{ {
0
(f
0
u
0
+ f
1
Du
0
+ · · · + f
p
D
p
u
0
)
=
{ {
0
(f
0
L + f
1
D + · · · + f
p
D
p
)u
0
=
Note u
0
= g D{
0
= D{ D{
0
= D({ {
0
) so that
{ {
p+1
= { {
0
(f
0
L + f
1
D + · · · + f
p
D
p
)D({ {
0
)
= (
L (f
0
L + f
1
D + · · · + f
p
D
p
))D({ {
0
)
= (
L (f
0
D + f
1
D + · · · + f
p
D
p+1
))({ {
0
)=
Thus
2(M({
p+1
) M({)) =
°
°
{
p+1
{
°
°
2
D
°
°
t
p+1
(D)({ {
0
)
°
°
2
D
(9.1.8)
where
t
p+1
(}) = 1 (f
0
} + f
1
}
2
+ · · · + f
p
}
p+1
)= One can make appropriate
choices of the polynomial
t
p+1
(}) and use some properties of eigenvalues and
Theorem 9.1.2 (CG Convergence Properties) Let
D be an q×q SPD matrix,
and consider
D{ = g=
1. The conjugate gradient method will obtain the solution within q iterations.
2. If g is a linear combination of n eigenvectors of D, then the conjugate
gradient method will obtain t he solution within n iterations.
3. If the set of all eigenvalues of
D has at most n distinct eigenvalues, then
the conjugate gradient method will obtain the solution within
n iterations.
4. Let
2
max
@
min
be the condition number of D> which is the ratio
of the largest and smallest eigenvalues of
D= Then the following error
estimate holds so that a condition number closer to one is desirable
k
{
p
{k
D
2
°
°
{
0
{
°
°
D
(
s
2
1
s
2
+ 1
)
p
= (9.1.9)
9.1.1 Exercises
1. Show the three prop erties if Krylov spaces in lines (9.1.1-9.1.3) are true.
2. Show equation (9.1.7) is true for
p = 2=
3. Give a mathematical induction proof of Theorem 9.1.1.
© 2004 by Chapman & Hall/CRC
matrix algebra to prove the following theorem, see [11, Chapter 2].
350 CHAPTER 9. KRYLOV METHODS FOR AX = D
9.2 Preconditioners
Consider the SPD matrix D and the problem D{ = g= For the preconditioned
conjugate gradient method a preconditioner is another SPD matrix P= The
matrix P is chosen so that the equivalent problem P
1
D{ = P
1
g can be
more easily solved. Note
P
1
D may not be SPD. However P
1
is also SPD
and must have a Cholesky factorization, which we write as
P
1
= V
W
V=
The system P
1
D{ = P
1
g can be rewritten in terms of a SPD matrix
b
Db{ =
b
g
V
W
VD{ = V
W
Vg
(VDV
W
)(V
W
{) = Vg where
b
D = VDV
W
, b{ = V
W
{ and
b
g = Vg=
b
D is SPD. It is similar to P
1
D, that
is,
V
W
(P
1
D)V
W
=
b
D so that the eigenvalues of P
1
D and
b
D are the same.
The idea is use the conjugate gradient method on
b
Db{ =
b
g and to choose P so
that the properties in Theorem 9.1.1 favor rapid convergence of the method.
At first glance it may appear to be very costly to do this. However the
following identities show that the application of the conjugate gradient method
to
b
Db{ =
b
g is relatively simple:
b
D = VDV
W
> b{ = V
W
{> bu = V(g D{) and bs = V
W
s=
Then
b
=
b
u
W
bu
bs
W
b
Dbs
=
(
Vu)
W
(Vu)
(V
W
s)
W
(VDV
W
) (V
W
s)
=
u
W
(V
W
Vu)
s
W
Ds
=
Similar calculations hold for
b
and bs. This eventually gives the following pre-
conditioned conjugate gradient algorithm, which does require an "easy" s olution
of
P} = u for each iterate.
Preconditioned Conjugate Gradient Method.
Let
{
0
be an initial guess
u
0
= g D{
0
solve P}
0
= u
0
and set s
0
= }
0
for m = 0, maxm
= (}
p
)
W
u
p
@(s
p
)
W
Ds
p
{
p+1
= {
p
+ s
p
u
p+1
= u
p
Ds
p
test for convergence
solve P}
p+1
= u
p+1
= (}
p+1
)
W
u
p+1
@(}
p
)
W
u
p
s
p+1
= }
p+1
+ s
p
endloop.
© 2004 by Chapman & Hall/CRC
9.2. PRECONDITIONERS 351
Examples of Preconditioners.
1. Block diagonal part of D= For the Poisson problem in two space variables
P =
5
9
9
9
7
E
0 0 0
0
E 0 0
0 0
.
.
.
.
.
.
0 0 · · ·
E
6
:
:
:
8
where
D =
5
9
9
9
7
E L
0 0
L E L 0
0
L
.
.
.
.
.
.
0 0 · · ·
E
6
:
:
:
8
and E =
5
9
9
9
7
4 1 0 0
1 4 1 0
0 1
.
.
.
.
.
.
0 0 · · · 4
6
:
:
:
8
=
2. Incomplete Cholesky factorization of
D= Let D = P + H =
b
J
b
J
W
+ H
where H is chosen so that either the smaller components in D are neglected
or some desirable structure of
P is attained. This can be done by defining a
subset
V {1> = = = > q} and overwriting the d
lm
when l> m 5 V
d
lm
= d
lm
d
ln
1
d
nn
d
nm
=
This common preconditioner and some variations are described in [6].
3. Incomplete domain decomposition. Let
D = P +H where P has the form
of an "arrow" associated with domain decompositions. For example, consider
a problem with two large blo cks separated by a smaller block, which is listed
last.
P =
5
7
D
11
0 D
13
0 D
22
D
23
D
31
D
32
D
33
6
8
where
D =
5
7
D
11
D
12
D
13
D
21
D
22
D
23
D
31
D
32
D
33
6
8
=
Related preconditioners are described in the paper by E. Chow and Y. Saad [5].
4. ADI sweeps. Form
P by doing one or more ADI sweeps. For example,
if
D = K + Y and one sweep is done in each direction, then for D = P Q as
defined by
(
L + K){
p+
1
2
= g + (L Y ){
p
(L + Y ){
p+1
= g + (L K){
p+
1
2
= g + (L K)((L + K)
1
(g + (L Y ){
p
)=
© 2004 by Chapman & Hall/CRC