White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI
Подождите немного. Документ загружается.
8.4. SOR AND P-REGULAR SPLITTINGS 331
Solve for
{ to obtain a fixed point problem
{ = P
1
Q{ + P
1
g=
The iterative method based on the splitting is
{
p+1
= P
1
Q{
p
+ P
1
g
= {
p
+ P
1
u({
p
)= (8.4.4)
Here
{
0
is some initial guess for the solution, and the solve step with the matrix
P is assumed to be relatively easy. The analysis of convergence can be either
in terms of some norm of
P
1
QD
an SPD matrix one can place conditions on the splitting so that the quadratic
function continues to decrease as the iteration advances. In particular, we will
show
M({
p+1
) = M({
p
)
1
2
(
{
p+1
{
p
)
W
(P
W
+ Q)({
p+1
{
p
)= (8.4.5)
So, if
P
W
+ Q is positive definite, then the sequence of real numbers M({
p
) is
decreasing. For more details on the following splitting consult J. Ortega [20].
Definition.
D = P Q is called a P-regular splitting if and only if P has
an inverse and P
W
+ Q is positive definite.
Note, if D is SPD, then D = D
W
= P
W
Q
W
and (P
W
+ Q )
W
= P + Q
W
=
P + P
W
D = P
W
+ Q= Thus, if D = P Q is P-regular and D is SPD,
then
P
W
+ Q is SPD.
Examples.
1. Jacobi splitting for
D = G (O + X) where P = G is the diagonal of D=
P
W
+ Q = G + (O + X) should be positive definite.
2. Gauss-Seidel splitting for
D = (G O) + X where P = G O is lower
triangular part of
D= P
W
+Q = (G O)
W
+X = G O
W
+X should be positive
definite. If
D is SPD, then O
W
= X and the G will have positive diagonal
components. In this case P
W
+ Q = G is SPD.
3. SOR splitting for 0
? $ ? 2 and D SPD with
D =
1
$
(G $O)
1
$
((1 $)G + $X)=
Here P =
1
$
(G $O) has an inverse because it is lower triangular with positive
diagonal components. This also gives a P-regular splitting because
P
W
+ Q = (
1
$
(G $O))
W
+
1
$
((1 $)G + $X)
= (
2
$
1)G=
© 2004 by Chapman & Hall/CRC
being less than one (see Section 2.5), or for
332 CHAPTER 8. CLASSICAL METHODS FOR AX = D
Theorem 8.4.2 (P-regular Splitting Convergence) If
D is SPD and D = P
Q
is a P-regular splitting, then the iteration in (8.4.4) will satisfy (8.4.5) and
will converge to the solution of
D{ = g.
Proof. First, the equality in (8.4.5) will be established by using (8.4.4) and
(8.4.3) with
| = {
p+1
= {
p
+ P
1
u({
p
) and { = {
p
M({
p+1
) = M({
p
) u({
p
)
W
P
1
u({
p
) +
1
2
(
P
1
u({
p
))
W
D(P
1
u({
p
))
=
M({
p
) u({
p
)
W
[P
1
1
2
P
W
DP
1
]u({
p
)
=
M({
p
) u({
p
)
W
P
W
[P
W
1
2
D]P
1
u({
p
)
= M({
p
) ({
p+1
{
p
)
W
[P
W
1
2
(
P Q)]({
p+1
{
p
)
=
M({
p
)
1
2
(
{
p+1
{
p
)
W
[2P
W
P + Q]({
p+1
{
p
)=
Since }
W
P
W
} = }
W
P}> (8.4.5) holds.
Second, we establish that
M({) is bounded from below. Use the inequality
in part 2 of Theorem 8.4.1,
{
W
D{ f
0
{
W
{, to write
M({) =
1
2
{
W
D{ {
W
g
1
2
f
0
{
W
{ {
W
g=
Next use the Cauchy inequality with D = L so that
¯
¯
{
W
g
¯
¯
¡
{
W
{
¢
1
2
¡
g
W
g
¢
1
2
,
and thus
M({)
1
2
f
0
{
W
{
¡
{
W
{
¢
1
2
¡
g
W
g
¢
1
2
1
2
f
0
[(
¡
{
W
{
¢
1
2
¡
g
W
g
¢
1
2
f
0
)
2
(
¡
g
W
g
¢
1
2
f
0
)
2
]
1
2
f
0
[0 (
¡
g
W
g
¢
1
2
f
0
)
2
]=
Third, note that M({
p
) is a decreasing sequence of real numbers that is
bounded from below. Since the real numbers are complete,
M({
p
) must con-
verge to some real number and M({
p
) M({
p+1
) =
1
2
({
p+1
{
p
)
W
(P
W
+
Q)({
p+1
{
p
) must converge to zero. Consider the norm associated with the
SPD matrix P
W
+ Q
°
°
{
p+1
{
p
°
°
2
P
W
+Q
= 2
¡
M({
p
) M({
p+1
)
¢
$ 0=
Thus, {
p+1
{
p
converges to the zero vector.
© 2004 by Chapman & Hall/CRC
8.5. SOR AND MPI 333
Fourth,
{
p+1
= {
p
+ P
1
u({
p
) so that P
1
u({
p
) converges to the zero
vector. Since
P is continuous, u({
p
) also converges to the zero vector. Since
D is SPD, there exists a solution of D{ = g, that is, u ({) = 0= Thus
u({
p
) u ({) = (g D{
p
) (g D{) = D({
p
{) $ 0=
Since D
1
is continuous, {
p
{ converges to the zero vector.
The next two sections will give additional examples of P-regular splittings
where the solve step P} = u({
p
) can be done in parallel. Such schemes can be
used as stand-alone algorithms or as preconditioners in the conjugate gradient
method.
8.4.1 Exercises
1. In the proof of part 3 in Theorem 8.4.1 prove the assertion about b
{
W
D|@|
W
D|=
2. Consider the following 2 × 2 SPD matrix and the indicated splitting
D =
E I
W
I F
¸
=
P
1
0
0
P
2
¸
Q
1
I
W
I Q
2
¸
=
(a). Find conditions on this splitting so that it will be P-regular.
(b). Apply this to SOR on the first and second blocks.
8.5 SOR and MPI
Consider a block 3 × 3 matrix with the following form
D =
5
7
D
11
D
12
D
13
D
21
D
22
D
23
D
31
D
32
D
33
6
8
=
One can think of this as a generalization of the Poisson problem with the first
two blocks of unknowns separated by a smaller i nterface block. When finite
di
erences or finite element methods are used, often the D
12
and D
21
are either
zero or sparse with small nonzero components. Consider splittings of the three
diagonal blocks
D
ll
= P
l
Q
l
=
Then one can associate a number of splittings with the large D matrix.
Examples.
1. In this example the
P in the splitting of D is block diagonal. The
inversion of
P has three independent computations
D =
5
7
P
1
0 0
0
P
2
0
0 0
P
3
6
8
5
7
Q
1
D
12
D
13
D
21
Q
2
D
23
D
31
D
32
Q
3
6
8
=
(8.5.1)
© 2004 by Chapman & Hall/CRC
334 CHAPTER 8. CLASSICAL METHODS FOR AX = D
2. The second example is slightly more complicated, and the
P is block
lower triangular. The inversion of
P has two independent computations and
some computations with just one processor.
D =
5
7
P
1
0 0
0 P
2
0
D
31
D
32
P
3
6
8
5
7
Q
1
D
12
D
13
D
21
Q
2
D
23
0 0 Q
3
6
8
=
(8.5.2)
3. This
P has the form from domain decomposition. The inversion of P
can be done by the Schur complement where E is the block diagonal of P
1
and
P
2
, and F = P
3
.
D =
5
7
P
1
0 D
13
0 P
2
D
23
D
31
D
32
P
3
6
8
5
7
Q
1
D
12
0
D
21
Q
2
0
0 0 Q
3
6
8
=
(8.5.3)
Theorem 8.5.1 (Block P-regular Splitting) Consider the iteration given by
(8.5.2). If
D is SPD, and
P
W
+ Q =
5
7
P
W
1
+ Q
1
D
12
0
D
21
P
W
2
+ Q
2
0
0 0
P
W
3
+ Q
3
6
8
is SPD, then the iteration will converge to the solution of D{ = g=
Proof. The proof is an easy application of Theorem 8.4.2. Since D is sym-
metric
D
W
lm
= D
ml
>
P
W
+ Q =
5
7
P
1
0 0
0 P
2
0
D
31
D
32
P
3
6
8
W
+
5
7
Q
1
D
12
D
13
D
21
Q
2
D
23
0 0 Q
3
6
8
=
5
7
P
W
1
0 D
W
31
0 P
W
2
D
W
32
0 0 P
W
3
6
8
+
5
7
Q
1
D
12
D
13
D
21
Q
2
D
23
0 0 Q
3
6
8
=
5
7
P
W
1
+ Q
1
D
12
0
D
21
P
W
2
+ Q
2
0
0 0
P
W
3
+ Q
3
6
8
=
If the D
12
= D
21
= 0> then P
W
+ Q will b e P-regular if and only if each
splitting
D
ll
= P
l
Q
l
is P-regular. A special case is the SOR algorithm
applied to the large blocks of unknowns, up dating the unknowns, and doing
SOR of the interface block, which is listed last. In this case P
l
=
1
$
(G
l
$O
l
)
so that
P
W
l
+ Q
l
= (
2
$
1)G
l
=
The following MPI code solves the Poisson problem where the spatial blocks
of nodes include horizontal large blocks with
q = 447 unknowns in the x direc-
tion and (
q s + 1) @s unknowns in the y direction. There are s = 2> 4> 8> 16
© 2004 by Chapman & Hall/CRC
8.5. SOR AND MPI 335
or 32 processors with (
q s + 1) @s = 223> 111> 55> 27 and 13 in unknowns in
the y direction, respectively. There are
s 1 smaller interface blocks with one
row each of
q unknowns so that the s + 1 block of the matrix D has q × (s 1)
unknowns=
The code is initialized in lines 6-30, and the SOR while loop is in lines 32-
120. SOR is done concurrently in lines 38-48 for the larger blocks, the results
are communicated to the adjacent processors in lines 50-83, then SOR is done
concurrently in lines 84-96 for the smaller interface blocks, and finally in lines
98-105 the results are communicated to adjacent blo cks. The communication
scheme used in lines 50-83 is similar to that used in Section 6.6 and illustrated
Here we have stored the
n
wk
big block and the top
interface block on processor
n 1 where 0 ? n ? s; the last big block is stored
on processor s 1= In lines 107-119 the global error is computed and broadcast
to all the processors. If it satisfies the convergence criteria, then the while loop
on each processor is exited. One can use the MPI subroutine mpi_allreduce()
code cgssormpi.f. The results in lines 123-133 are gathered onto processor 0
and partial results are printed.
MPI/Fortran Code sorddmpi.f
1. program sor
2.!
3.! Solve Poisson equation via SOR.
4.! Uses domain decomposition to attain parallel computation.
5.!
6. implicit none
7. include ’mpif.h’
8. real ,dimension (449,449)::u,uold
9. real ,dimension (1:32)::errora
10. real :: w, h, eps,pi,error,utemp,to,t1,timef
11. integer :: n,maxk,maxit,it,k,i,j,jm,jp
12. integer :: my_rank,p,source,dest,tag,ierr,loc_n
13. integer :: status(mpi_status_size),bn,en,sbn
14. n = 447
15. w = 1.99
16. h = 1.0/(n+1)
17. u = 0.
18. errora(:) = 0.0
19. error = 1.
20. uold = 0.
21. call mpi_init(ierr)
22. call mpi_comm_rank(mpi_comm_world,my_rank,ierr)
23. call mpi_comm_size(mpi_comm_world,p,ierr)
24. if (my_rank.eq.0) then
25. to = timef()
© 2004 by Chapman & Hall/CRC
in Figures 6.6.1 and 6.6.2.
to combine the gather and broadcast operations, see Section 9.3 and the MPI
336 CHAPTER 8. CLASSICAL METHODS FOR AX = D
26. end if
27. pi = 3.141592654
28. maxit = 2000
29. eps = .001
30. it = 0
31.! Begin the while loop for the parallel SOR iterations.
32. do while ((it.lt.maxit).and.(error.gt.eps))
33. it = it + 1
34. loc_n = (n-p+1)/p
35. bn = 2+(my_rank)*(loc_n+1)
36. en = bn + loc_n -1
37.! Do SOR for big blocks.
38. do j=bn,en
39. do i =2,n+1
40. utemp = (1000.*sin((i-1)*h*pi)*sin((j-1)*h*pi)*h*h&
41. + u(i-1,j) + u(i,j-1) &
42. + u(i+1,j) + u(i,j+1))*.25
43. u(i,j) = (1. -w)*u(i,j) + w*utemp
44. end do
45. end do
46. errora(my_rank+1) = maxval(abs(u(2:n+1,bn:en)-&
47. uold(2:n+1,bn:en)))
48. uold(2:n+1,bn:en) = u(2:n+1,bn:en)
49.! Communicate computations to adjacent blocks.
50. if (my_rank.eq.0) then
51. call mpi_recv(u(1,en+2),(n+2),mpi_real,my_rank+1,50,&
52. mpi_comm_world,status,ierr)
53. call mpi_send(u(1,en+1),(n+2),mpi_real,my_rank+1,50,&
54. mpi_comm_world,ierr)
55. end if
56. if ((my_rank.gt.0).and.(my_rank.lt.p-1)&
57. .and.(mod(my_rank,2).eq.1)) then
58. call mpi_send(u(1,en+1),(n+2),mpi_real,my_rank+1,50,&
59. mpi_comm_world,ierr)
60. call mpi_recv(u(1,en+2),(n+2),mpi_real,my_rank+1,50,&
61. mpi_comm_world,status,ierr)
62. call mpi_send(u(1,bn),(n+2),mpi_real,my_rank-1,50,&
63. mpi_comm_world,ierr)
64. call mpi_recv(u(1,bn-1),(n+2),mpi_real,my_rank-1,50,&
65. mpi_comm_world,status,ierr)
66. end if
67. if ((my_rank.gt.0).and.(my_rank.lt.p-1)&
68. .and.(mod(my_rank,2).eq.0)) then
69. call mpi_recv(u(1,bn-1),(n+2),mpi_real,my_rank-1,50,&
70. mpi_comm_world,status,ierr)
© 2004 by Chapman & Hall/CRC
8.5. SOR AND MPI 337
71. call mpi_send(u(1,bn),(n+2),mpi_real,my_rank-1,50,&
72. mpi_comm_world,ierr)
73. call mpi_recv(u(1,en+2),(n+2),mpi_real,my_rank+1,50,&
74. mpi_comm_world,status,ierr)
75. call mpi_send(u(1,en+1),(n+2),mpi_real,my_rank+1,50,&
76. mpi_comm_world,ierr)
77. end if
78. if (my_rank.eq.p-1) then
79. call mpi_send(u(1,bn),(n+2),mpi_real,my_rank-1,50,&
80. mpi_comm_world,ierr)
81. call mpi_recv(u(1,bn-1),(n+2),mpi_real,my_rank-1,50,&
82. mpi_comm_world,status,ierr)
83. end if
84. if (my_rank.lt.p-1) then
85. j = en +1
86.! Do SOR for smaller interface blocks.
87. do i=2,n+1
88. utemp = (1000.*sin((i-1)*h*pi)*sin((j-1)*h*pi)*h*h&
89. + u(i-1,j) + u(i,j-1)&
90. + u(i+1,j) + u(i,j+1))*.25
91. u(i,j) = (1. -w)*u(i,j) + w*utemp
92. end do
93. errora(my_rank+1) = max1(errora(my_rank+1),&
94. maxval(abs(u(2:n+1,j)-uold(2:n+1,j))))
95. uold(2:n+1,j) = u(2:n+1,j)
96. endif
97.! Communicate computations to adjacent blocks.
98. if (my_rank.lt.p-1) then
99. call mpi_send(u(1,en+1),(n+2),mpi_real,my_rank+1,50,&
100. mpi_comm_world,ierr)
101. end if
102. if (my_rank.gt.0) then
103. call mpi_recv(u(1,bn-1),(n+2),mpi_real,my_rank-1,50,&
104. mpi_comm_world,status,ierr)
105. end if
106.! Gather local errors to pro cessor 0.
107. call mpi_gather(errora(my_rank+1),1,mpi_real,&
108. errora,1,mpi_real,0,&
109. mpi_comm_world,ierr)
110. call mpi_barrier(mpi_comm_world,ierr)
111.! On processor 0 compute the maximum of the local errors.
112. if (my_rank.eq.0) then
113. error = maxval(errora(1:p))
114. end if
115.! Send this global error to all processors so that
© 2004 by Chapman & Hall/CRC
338 CHAPTER 8. CLASSICAL METHODS FOR AX = D
Table 8.5.1: MPI Times for sorddmpi.f
p time iteration
2 17.84 673
4 07.93 572
8 03.94 512
16 02.18 483
32 01.47 483
116.! they will exit the while loop when the gl obal error
117.! test is satisfied.
118. call mpi_bcast(error,1,mpi_real,0,&
119. mpi_comm_world,ierr)
120. end do
121.! End of the while loop.
122.! Gather the computations to processor 0
123. call mpi_gather(u(1,2+my_rank*(loc_n+1)),&
124. (n+2)*(loc_n+1),mpi_real,&
125. u,(n+2)*(loc_n+1),mpi_real,0,&
126. mpi_comm_world,ierr)
127. if (my_rank.eq.0) then
128. t1 = timef()
129. print*, ’sor iterations = ’,it
130. print*, ’time = ’, t1
131. print*, ’error = ’, error
132. print*, ’center value of solution = ’, u(225,225)
133. end if
134. call mpi_finalize(ierr)
135. end
The computations in Table 8.5.1 use 2, 4, 8, 16 and 32 processors. Since the
orderings of the unknowns are di
erent, the SOR iterations required for conver-
gence vary; in this case they decrease as the number of processors increase. Also,
as the number of processors increase, the number of interface blocks increase
as does the amount of communication. The execution times (sec.) reflect rea-
sonable speedups given the decreasing iterations for convergence and increasing
communication.
8.5.1 Exercises
1. Establish analogues of Theorem 8.5.1 for the splittings in lines (8.5.1)
and (8.5.3).
2. In Theorem 8.5.1 use the Schur complement to give conditions so that
P
W
+ Q will be positive definite.
© 2004 by Chapman & Hall/CRC
8.6. PARALLEL ADI SCHEMES 339
3.
criteria.
4. Experiment with sorddmpi.f by varying the number of unknowns and
SOR parameter. Note the e
ects of using dierent numbers of processors.
8.6 Parallel ADI Schemes
Alternating direction implicit (ADI) iterative methods can be used to approxi-
mate the solution to a two variable Poisson problem by a sequence of ordinary
dierential equations solved in the x or y directions
x
{{
x
||
= i
x
{{
= i + x
||
x
||
= i + x
{{
=
One may attempt to approximate this by the following scheme
for m = 0, maxm
for each y solve
x
p+
1
2
{{
= i + x
p
||
for each x solve
x
p+1
||
= i + x
p+
1
2
{{
test for convergence
endloop.
The discrete version of this uses the finite di
erence method to discretize
the ordinary di
erential equations
x
p+
1
2
l+1>m
2x
p+
1
2
l>m
+ x
p+
1
2
l1>m
{
2
= i
l>m
+
x
p
l>m
+1
2x
p
l>m
+ x
p
l>m
1
|
2
(8.6.1)
x
p+1
l>m+1
2x
p+1
l>m
+ x
p+1
l>m1
|
2
= i
l>m
+
x
p+
1
2
l+1>m
2x
p+
1
2
l>m
+ x
p+
1
2
l1>m
{
2
= (8.6.2)
The discrete version of the above scheme is
for m = 0, maxm
for each
m solve
tridiagonal problem (8.6.1) for
x
p+
1
2
l>m
for each l solve
tridiagonal problem (8.6.2) for
x
p+1
l>m
test for convergence
endloop.
The solve steps may be done by the tridiagonal algorithm. The
l or m solves are
independent, that is, the solves for
x
p+
1
2
l>m
and x
p+
1
2
l>
b
m
can be done concurrently
© 2004 by Chapman & Hall/CRC
Verify the computations in Table 8.5.1. Exp eriment with the convergence
340 CHAPTER 8. CLASSICAL METHODS FOR AX = D
for
m 6=
b
m, and also the solves for x
p+1
l>m
and x
p+1
b
l>m
can be done concurrently for
l 6=
b
l.
The matrix form of this scheme requires the coe
!cient matrix to be written
as
D = K + Y where K and Y reflect the discretized ordinary dierential
equations in the x and y directions, respectively. For example, suppose { =
| = k and there are q unknowns in each direction of a square domain. The
D will be a block q × q matrix and using the classical order of bottom grid row
first and moving from left to right in the grid rows
D = K + Y
5
7
E L ]
L E L
] L E
6
8
=
5
7
F ] ]
] F ]
] ] F
6
8
+
5
7
2L L ]
L
2L L
] L
2L
6
8
where L is a q × q identity matrix, ] is a q ×
q zero matrix and F is a q × q
tridiagonal matrix, for example for q = 3
F =
5
7
2 1 0
1 2 1
0
1 2
6
8
=
The ADI algorithm for the solution of D{ = g is based on two splittings of
D
D = (L + K) (L Y )
D = (L + Y ) (L K)=
The positive constant is chosen so that L + K and L + Y have inverses and
to accelerate convergence. More generally,
L is replaced by a diagonal matrix
with positive diagonal components chosen to obtain an optimal convergence
rate.
ADI Algorithm for D{ = g with D = K + Y=
for m = 0, maxm
solve (
L + K){
p+
1
2
= g + (L Y ){
p
solve (L + Y ){
p+1
= g + (L K){
p+
1
2
test for convergence
endloop.
This algorithm may be written in terms of a single splitting.
{
p+1
= (L + Y )
1
[g + (L K){
p+
1
2
]
= (
L + Y )
1
[g + (L K)(L + K)
1
(g + (L Y ){
p
)]
= (
L + Y )
1
[g + (L K)(L + K)
1
g] +
(
L + Y )
1
(L K)(L + K)
1
(L Y ){
p
= P
1
g + P
1
Q{
p
© 2004 by Chapman & Hall/CRC