White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI
Подождите немного. Документ загружается.
352 CHAPTER 9. KRYLOV METHODS FOR AX = D
Solve for
{
p+1
= P
1
g + P
1
Q{
p
where
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
so that
P = (L + K)
1
2
(L + Y )=
The parallelism in this approach comes from independent tridiagonal solves in
each direction. Another approach is to partition each tridiagonal solve as is
indicated in [15].
5. SSOR. Form
P by doing one forward SOR and then one backward SOR
sweep. Let
D = G O X where O
W
= X and 0 ? $ ? 2=
1
$
(G $O){
p+
1
2
= g +
1
$
((1 $)G + $X){
p
1
$
(G $X ){
p+1
= g +
1
$
((1 $)G + $O){
p+
1
2
= g +
1
$
((1 $)G + $O)((
1
$
(G $O))
1
(g +
1
$
((1 $)G + $X){
p
)=
Solve for {
p+1
= P
1
g + P
1
Q{
p
where
P
1
= (
1
$
(G $X ))
1
(L +
1
$
((1 $)G + $O)(
1
$
(G $O))
1
= (
1
$
(G $X ))
1
((
1
$
(G $O)) +
1
$
((1 $)G + $O)(
1
$
(G $O))
1
= (
1
$
(G $X ))
1
2 $
$
G
(
1
$
(G $O))
1
so that
P =
1
$
(G $O)(
2
$
$
G
)
1
1
$
(G $X )=
6. Additive Schwarz. As motivation consider the Poisson problem and di-
vide the space grid into subsets of q
l
unknowns so that q
l
q= Let E
l
be
associated with a splitting of a restriction of the coe
!cient matrix to subset l
of unknowns, and let U
l
be the restriction operator so that
D : R
q
$ R
q
E
l
: R
q
l
$ R
q
l
U
l
: R
q
$ R
q
l
U
W
l
: R
q
l
$ R
q
=
Define
c
P
l
= U
W
l
E
1
l
U
l
: R
q
$ R
q
= Although these matrices are not invertible,
one may be able to associate a splitting with the summation
c
P
l
. Often a
© 2004 by Chapman & Hall/CRC
9.2. PRECONDITIONERS 353
coarse mesh is associated with the problem. In this case let
D
0
b e q
0
× q
0
and
let
U
W
0
: R
q
0
$ R
q
be an extension operator from the coarse to fine mesh.
The following may accelerate the convergence
P
1
U
W
0
D
1
0
U
0
+
c
P
l
= U
W
0
D
1
0
U
0
+ U
W
l
E
1
l
U
l
=
A common example is to apply SSOR to the subsets of unknowns with zero
boundary conditions. A nice survey of Schwarz methods can be found in the
pap er written by Xiao-Chuan Cai in the first chapter of [12].
7. Least s quares approximations for
D
c
P = L= This is equivalent to q least
squares problems
D bp
m
= h
m
=
For example, if the column vectors bp
m
were to have nonzero components in
rows
l = m 1> m and m + 1> then this becomes a least squares problem with q
equations and three unknowns ep
m
= [p
m1
> p
m
> p
m+1
]
W
where D is restricted
to columns
m 1> m and m + 1
D(1 : q> m 1 : m + 1) ep
m
= h
m
=
The preconditioner is formed by collecting the column vectors bp
m
P
1
= [
b
p
1
bp
2
· · · bp
q
]=
Additional information on this approach can be found in the pap er by M. J.
Grote and T. Huckle [9].
The following M
ATLA B code is a slight variation of precg.m that was de-
scribed in Section 3.6. Here the SSOR and the block diagonal preconditioners
are used. The choice of preconditioners is made in lines 28-33. The number of
iterates required for convergence was 19, 55 and 73 for SSOR, block diagonal
and no preconditioning, respectively. In the ssorpc.m preconditioner function
the forward solve is done in lines 3-7, and the backward solve is done in lines
9-13. In the bdiagpc.m preconditioner function the diagonal blocks are all the
same and are defined in lines 3-12 every time the function is evaluated. The
solves for each block are done in lines 13-15.
MATLAB Code pccg.m with ssorpc.m and bdiagpc.m
1. % Solves -u_xx - u_yy = 200+200sin(pi x)sin(pi y).
2. % Uses PCG with SSOR or block diagonal preconditioner.
3. % Uses 2D arrays for the column vectors.
4. % Does not explicity store the matrix.
5. clear;
6. w = 1.6;
7. n = 65;
8. h = 1./n;
9. u(1:n+1,1:n+1)= 0.0;
© 2004 by Chapman & Hall/CRC
354 CHAPTER 9. KRYLOV METHODS FOR AX = D
10. r(1:n+1,1:n+1)= 0.0;
11. rhat(1:n+1,1:n+1) = 0.0;
12. % Define right side of PDE
13. for j= 2:n
14. for i = 2:n
15. r(i,j)= h*h*(200+200*sin(pi*(i-1)*h)*sin(pi*(j-1)*h));
16. end
17. end
18. errtol = .0001*sum(sum(r(2:n,2:n).*r(2:n,2:n)))^.5;
19. p(1:n+1,1:n+1)= 0.0;
20. q(1:n+1,1:n+1)= 0.0;
21. err = 1.0;
22. m = 0;
23. rho = 0.0;
24. % Begin PCG iterations
25. while ((err
A errtol)&(m ? 200))
26. m = m+1;
27. oldrho = rho;
28. % Execute SSOR preconditioner
29. rhat = ssorpc(n,n,1,1,1,1,4,.25,w,r,rhat);
30. % Execute block diagonal preconditioner
31. % rhat = bdiagpc(n,n,1,1,1,1,4,.25,w,r,rhat);
32. % Use the following line for no preconditioner
33. % rhat = r;
34. % Find conjugate direction
35. rho = sum(sum(r(2:n,2:n).*rhat(2:n,2:n)));
36. if (m==1)
37. p = rhat;
38. else
39. p = rhat + (rho/oldrho)*p;
40. end
41. % Use the following line for steepest descent method
42. % p=r;
43. % Executes the matrix product q = Ap without storage of A
44. for j= 2:n
45. for i = 2:n
46. q(i,j)=4.*p(i,j)-p(i-1,j)-p(i,j-1)-p(i+1,j)-p(i,j+1);
47. end
48. end
49. % Executes the steepest descent segment
50. alpha = rho/sum(sum(p.*q));
51. u = u + alpha*p;
52. r = r - alpha*q;
53. % Test for convergence via the infinity norm of the residual
54. err = max(max(abs(r(2:n,2:n))));
© 2004 by Chapman & Hall/CRC
9.2. PRECONDITIONERS 355
55. reserr(m) = err;
56. end
57. m
58. semilogy(reserr)
1. function r = ssorpc(nx,ny,ae,aw,as,an,ac,rac,w,d,r)
2. % This preconditioner is SSOR.
3. for j= 2:ny
4. for i = 2:nx
5. r(i,j) = w*(d(i,j) + aw*r(i-1,j) + as*r(i,j-1))*rac;
6. end
7. end
8. r(2:nx,2:ny) = ((2.-w)/w)*ac*r(2:nx,2:ny);
9. for j= ny:-1:2
10. for i = nx:-1:2
11. r(i,j) = w*(r(i,j)+ae*r(i+1,j)+an*r(i,j+1))*rac;
12. end
13. end
1. function r = bdiagp c(nx,ny,ae,aw,as,an,ac,rac,w,d,r)
2. % This preconditioner is block diagonal.
3. Adiag = zeros(nx-1);
4. for i = 1:nx-1
5. Adiag(i,i) = ac;
6. if i
A1
7. Adiag(i,i-1) = -aw;
8. end
9. if i
?nx-1
10. Adiag(i,i+1) = -ae;
11. end
12. end
13. for j = 2:ny
14. r(2:nx,j) = Adiag\d(2:nx,j);
15. end
9.2.1 Exercises
1. In the derivation of the preconditioned conjugate gradient method do the
calculations for
b
and bs= Complete the derivation of the PCG algorithm.
2. Verify the calculations for the MATLA B code pccg.m. Experiment with
some variations of the SSOR and block preconditioners.
3. Experiment with the incomplete Cholesky preconditioner.
4. Experiment with the incomplete domain decomp osition preconditioner.
5. Experiment with the ADI preconditioner.
6. Experiment with the additive Schwarz preconditioner.
© 2004 by Chapman & Hall/CRC
356 CHAPTER 9. KRYLOV METHODS FOR AX = D
7. Experiment with the least squares preconditioner. You may want to
review least squares as in Section 9.4 and see the interesting paper by Grote
9.3 PCG and MPI
This section contains the MPI Fortran code of the preconditioned conjugate
gradient algorithm for the solution of a Poisson problem. The preconditioner is
an implementation of the additive Schwarz preconditioner with no coarse mesh
acceleration. It uses SSOR on large blocks of unknowns by partitioning the
second index with zero boundary conditions on the grid boundaries. So, this
could be viewed as a block diagonal preconditioner where the diagonal blocks
of the coe
!cient matrix are split by the SSOR splitting. Since the number of
blocks are associated with the number of processors, the preconditioner really
does change with the number of processors.
The global initialization is done in lines 1-36, and an initial guess is the
zero vector. In lines 37-42 MPI is started and the second index is partitioned
according to the number of processors. The conjugate gradient loop is executed
in lines 48-120, and the partial outputs are given by lines 121-138.
The conjugate gradient loop has substages, which are done in parallel. The
preconditioner is done on each block in lines 50-62. The local dot products for
are computed, then mpi_allreduce() is used to total the local dot products
and to broadcast the result to all the processors. The local parts of the updated
search direction are computed in lines 67-71. In order to do the sparse matrix
product
Ds, the top and bottom grid rows of s are communicated in lines 72-
107 to the adjacent processors. This communication scheme is similar to that
contain the local computation of
Ds. In lines 111-114 the local dot products
for the computation of
in the steepest descent direction computation are
computed, and then mpi_allreduce() is used to total the local dot products and
to broadcast the result to all the processors. Lines 114 and 115 contain the
updated local parts of the approximated solution and residual. Lines 117-118
contain the local computation of the residual error, and then mpi_allreduce() is
used to total the errors and to broadcast the result to all the processors. Once
the error criteria have been satisfied for all processors, the conjugate gradient
loop will be exited.
MPI/Fortran Code cgssormpi.f
1. program cgssor
2.! This code approximates the solution of
3.! -u_xx - u _yy = f
4.! PCG is used with a SSOR verson of the
5.! Schwarz additive preconditioner.
6.! The sparse matrix product, dot products and updates
7.! are also done in parallel.
© 2004 by Chapman & Hall/CRC
used in Section 6.6 and is illustrated in Figures 6.6.1 and 6.6.2.Lines 108-109
and Huckle [9].
9.3. PCG AND MPI 357
8. implicit none
9. include ’mpif.h’
10. real,dimension(0:1025,0:1025):: u,p,q,r,rhat
11. real,dimension (0:1025) :: x,y
12. real :: oldrho,ap, rho,alpha,error,dx2,w,t0,timef,tend
13. real :: loc_rho,loc_ap,loc_error
14. integer :: i,j,n,m
15. integer :: my_rank,proc,source,dest,tag,ierr,loc_n
16. integer :: status(mpi_status_size),bn,en
17. integer :: maxit,sbn
18. w = 1.8
19. u = 0.0
20. n = 1025
21. maxit = 200
22. dx2 = 1./(n*n)
23. do i=0,n
24. x(i) = float(i)/n
25. y(i) = x(i)
26. end do
27. r = 0.0
28. rhat = 0.0
29. q = 0.0
30. p = 0.0
31. do j = 1,n-1
32. r(1:n-1,j)=200.0*dx2*(1+sin(3.14*x(1:n-1))*sin(3.14*y(j)))
33. end do
34. error = 1.
35. m = 0
36. rho = 0.0
37. call mpi_init(ierr)
38. call mpi_comm_rank(mpi_comm_world,my_rank,ierr)
39. call mpi_comm_size(mpi_comm_world,proc,ierr)
40. loc_n = (n-1)/proc
41. bn = 1+(my_rank)*loc_n
42. en = bn + loc_n -1
43. call mpi_barrier(mpi_comm_world,ierr)
44. if (my_rank.eq.0) then
45. t0 = timef()
46. end if
47. do while ((error
A.0001).and.(m?maxit))
48. m = m+1
49. oldrho = rho
50. ! Execute Schwarz additive SSOR preconditioner.
51. ! This preconditioner changes with the number of processors!
52 do j= bn,en
© 2004 by Chapman & Hall/CRC
358 CHAPTER 9. KRYLOV METHODS FOR AX = D
53. do i = 1,n-1
54. rhat(i,j) = w*(r(i,j)+rhat(i-1,j)+rhat(i,j-1))*.25
55. end do
56. end do
57. rhat(1:n-1,bn:en) = ((2.-w)/w)*4.*rhat(1:n-1,bn:en)
58. do j= en,bn,-1
59. do i = n-1,1,-1
60. rhat(i,j) = w*(rhat(i,j)+rhat(i+1,j)+rhat(i,j+1))*.25
61. end do
62. end do
63. ! rhat = r
64. ! Find conjugate direction.
65. loc_rho = sum(r(1:n-1,bn:en)*rhat(1:n-1,bn:en))
66. call mpi_allreduce(loc_rho,rho,1,mpi_real,mpi_sum,&
mpi_comm_world,ierr)
67. if (m.eq.1) then
68. p(1:n-1,bn:en) = rhat(1:n-1,bn:en)
69. else
70. p(1:n-1,bn:en) = rhat(1:n-1,bn:en)&
+ (rho/oldrho)*p(1:n-1,bn:en)
71. endif
72. ! Execute matrix product q = Ap.
73. ! First, exchange information between processors.
74. if (my_rank.eq.0) then
75. call mpi_recv(p(0,en+1),(n+1),mpi_real,my_rank+1,50,&
76. mpi_comm_world,status,ierr)
77. call mpi_send(p(0,en),(n+1),mpi_real,my_rank+1,50,&
78. mpi_comm_world,ierr)
79. end if
80. if ((my_rank.gt.0).and.(my_rank.lt.proc-1)&
81. .and.(mod(my_rank,2).eq.1)) then
82. call mpi_send(p(0,en),(n+1),mpi_real,my_rank+1,50,&
83. mpi_comm_world,ierr)
84. call mpi_recv(p(0,en+1),(n+1),mpi_real,my_rank+1,50,&
85. mpi_comm_world,status,ierr)
86. call mpi_send(p(0,bn),(n+1),mpi_real,my_rank-1,50,&
87. mpi_comm_world,ierr)
88. call mpi_recv(p(0,bn-1),(n+1),mpi_real,my_rank-1,50,&
89. mpi_comm_world,status,ierr)
90. end if
91. if ((my_rank.gt.0).and.(my_rank.lt.proc-1)&
92. .and.(mod(my_rank,2).eq.0)) then
93. call mpi_recv(p(0,bn-1),(n+1),mpi_real,my_rank-1,50,&
94. mpi_comm_world,status,ierr)
95. call mpi_send(p(0,bn),(n+1),mpi_real,my_rank-1,50,&
© 2004 by Chapman & Hall/CRC
9.3. PCG AND MPI 359
96. mpi_comm_world,ierr)
97. call mpi_recv(p(0,en+1),(n+1),mpi_real,my_rank+1,50,&
98. mpi_comm_world,status,ierr)
99. call mpi_send(p(0,en),(n+1),mpi_real,my_rank+1,50,&
100. mpi_comm_world,ierr)
101. end if
102. if (my_rank.eq.proc-1) then
103. call mpi_send(p(0,bn),(n+1),mpi_real,my_rank-1,50,&
104. mpi_comm_world,ierr)
105. call mpi_recv(p(0,bn-1),(n+1),mpi_real,my_rank-1,50,&
106. mpi_comm_world,status,ierr)
107. end if
108. q(1:n-1,bn:en)=4.0*p(1:n-1,bn:en)-p(0:n-2,bn:en)-p(2:n,bn:en)&
109. - p(1:n-1,bn-1:en-1) - p(1:n-1,bn+1:en+1)
110. ! Find steepest descent.
111. loc_ap = sum(p(1:n-1,bn:en)*q(1:n-1,bn:en))
112. call mpi_allreduce(loc_ap,ap,1,mpi_real,mpi_sum,&
113. mpi_comm_world,ierr)
114. alpha = rho/ap
115. u(1:n-1,bn:en) = u(1:n-1,bn:en) + alpha*p(1:n-1,bn:en)
116. r(1:n-1,bn:en) = r(1:n-1,bn:en) - alpha*q(1:n-1,bn:en)
117. loc_error = maxval(abs(r(1:n-1,bn:en)))
118. call mpi_allreduce(loc_error,error,1,mpi_real, mpi_sum,&
119. mpi_comm_world,ierr)
120. end do
121. ! Send local solutions to processor zero.
122. if (my_rank.eq.0) then
123. do source = 1,proc-1
124. sbn = 1+(source)*loc_n
125. call mpi_recv(u(0,sbn),(n+1)*loc_n,mpi_real,source,50,&
126. mpi_comm_world,status,ierr)
127. end do
128. else
129. call mpi_send(u(0,bn),(n+1)*loc_n,mpi_real,0,50,&
130. mpi_comm_world,ierr)
131. end if
132. if (my_rank.eq.0) then
133. tend = timef()
134. print*, ’time =’, tend
135. print*, ’time per iteration = ’, tend/m
136. print*, m,error, u(512 ,512)
137. print*, ’w = ’,w
138. end if
139. call mpi_finalize(ierr)
140. end program
© 2004 by Chapman & Hall/CRC
360 CHAPTER 9. KRYLOV METHODS FOR AX = D
Table 9.3.1: MPI Times for cgssormpi.f
p time iteration
2 35.8 247
4 16.8 260
8 07.9 248
16 03.7 213
32 03.0 287
The Table 9.3.1 contains computations for q = 1025 and using z = 1=8= The
computation times are in seconds, and note the number of iterations vary with
the number of processors.
9.3.1 Exercises
1. Verify the computations in Table 9.3.1. Experiment with the convergence
criteria.
2. Experiment with variations on the SSOR preconditioner and include dif-
ferent
q and $=
3. Experiment with variations of the SSOR preconditioner to include the
use of a coarse mesh in the additive Schwarz preconditioner.
4. Use an ADI preconditioner in place of the SSOR preconditioner.
9.4 Least Squares
Consider an algebraic system where there are more equations than unknowns.
This will be a subproblem in the next two sections where the unknowns will be
the coe
!cients of the Krylov vectors. Let D be q × p where q A p= In this
case it may not be possible to find
{ such that D{ = g> that is, the residual
vector
u({) = g D{ may never be the zero vector. The next best alternative
is to find { so that in some way the residual vector is as small as possible.
Definition. Let
U({) u({)
W
u({) where D is q × p> u({) = g D{ and { is
p × 1= The least squares solution of D{ = g is
U({) = min
|
U(|)=
The following identity is important in finding a least squares solution
U(|) = (g D|)
W
(g D|)
=
g
W
g 2(D|)
W
g + (D|)
W
D|
= g
W
g + 2[1@2 |
W
(D
W
D)| |
W
(D
W
g)]= (9.4.1)
If
D
W
D is SPD, then by Theorem 8.4.1 the second term in (9.4.1) will be a
minimum if and only if
D
W
D{ = D
W
g= (9.4.2)
© 2004 by Chapman & Hall/CRC
9.4. LEAST SQUARES 361
This system is called the normal equations.
Theorem 9.4.1 (Normal Equations) If D has full column rank (D{ = 0 im-
plies
{ = 0), then the least squares solution is characterized by the solution of
the normal equations (9.4.2).
Proof. Clearly
D
W
D is symmetric. Note {
W
(D
W
D){ = (D{)
W
(D{) = 0 if
and only if D{ = 0= The full column rank assumption implies { = 0 so that
{
W
(D
W
D){ A 0 if { 6= 0= Thus D
W
D is SPD. Apply the first part of Theorem
8.4.1 to the second term in (9.4.1). Since the first term in (9.4.1) is constant
with respect to |> U(|) will be minimized if and only if the normal equations
(9.4.2) are satisfied.
Example 1. Consider the 3 × 2 algebraic system
5
7
1 1
1 2
1 3
6
8
{
1
{
2
¸
=
5
7
4
7
8
6
8
=
This could have evolved from the linear curve | = pw + f fit to the data
(w
l
> |
l
) = (1> 4)> (2> 7) and (3> 8) where {
1
= f and {
2
= p= The matrix has full
column rank and the normal equations are
3 6
6 14
¸
{
1
{
2
¸
=
19
42
¸
=
The solution is {
1
= f = 7@3 and {
2
= p = 2.
The normal equations are often ill-conditioned and prone to significant ac-
cumulation of roundo
errors. A good alternative is to use a QR factorization
of A.
Definition. Let
D be q × p= Factor D = TU where T is q × p such that
T
W
T = L> and U is p×p is upper triangular. This is called a QR factorization
of D=
Theorem 9.4.2 (QR Factorization) If D = TU and has full column rank,
then the solution of the normal equations is given by the solution of
U{ = T
W
g=
Proof. The normal equations become
(
TU)
W
(TU){ = (TU)
W
g
U
W
(T
W
T)U{ = U
W
T
W
g
U
W
U{ = U
W
T
W
g=
Because D is assumed to have full column rank, U must have an inverse. Thus
we only need to solve
U{ = T
W
g=
There are a number of ways to find the QR factorization of the matrix.
The modified Gram-Schmidt method is often used when the matrix has mostly
© 2004 by Chapman & Hall/CRC