White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI
Подождите немного. Документ загружается.
310 CHAPTER 7. MESSAGE PASSING INTERFACE
47. !
48. subroutine Fox(n,grid,local_A,local_B,local_C,n_bar)
49. integer, intent(in) :: n, n_bar
50. TYPE(GRID_INFO_TYPE), intent(in) :: grid
51. real, intent(in) , dimension(:,:) :: local_A, local_B
52. real, intent(out), dimension (:,:):: local_C
53. real, dimension(1:n_bar,1:n_bar) :: temp_A
54. integer:: step, source, dest, request,i,j
55. integer:: status(MPI_STATUS_SIZE), bcast_root
56. temp_A = 0.0
57. local_C = 0.0
58. source = mod( (grid%my_row + 1), grid%q )
59. dest = mod( (grid%my_row - 1 + grid%q), grid%q )
60. do step = 0, grid%q -1
61. bcast_root = mod( (grid%my_row + step), grid%q )
62. if (bcast_root == grid%my_col) then
63. call mpi_bcast(local_A, n_bar*n_bar, mpi_real,&
64. bcast_root, grid%row_comm, ierr)
65. ! print*, grid%my_row, grid%my_col, ’local_A = ’,local_A
66. call sgemm(’N’,’N’,n_bar,n_bar,n_bar,1.0,&
67. local_A,n_bar,local_B,n_bar,1.0,local_C,n_bar)
68. ! do j = 1,n_bar
69. ! do k = 1,n_bar
70. ! do i = 1,n_bar
71. ! local_C(i,j)=local_C(i,j) + local_A(i,k)*&
72. ! local_B(k,j)
73. ! end do
74. ! end do
75. ! end do
76. else
77. ! Store local_A from bcast_root in temp_A so as
78. ! not to overwrite local_A in destination.
79. call mpi_bcast(temp_A, n_bar*n_bar, mpi_real,&
80. bcast_root, grid%row_comm, ierr)
81. call sgemm(’N’,’N’,n_bar,n_bar,n_bar,1.0,&
82. temp_A,n_bar,local_B,n_bar,1.0,local_C,n_bar)
83. ! do j = 1,n_bar
84. ! do k = 1,n_bar
85. ! do i = 1,n_bar
86. ! local_C(i,j)=local_C(i,j) + temp_A(i,k)*&
87. ! local_B(k,j)
88. ! enddo
89. ! enddo
90. ! enddo
91. endif
© 2004 by Chapman & Hall/CRC
7.6. FOX ALGORITHM FOR AB 311
Table 7.6.1: Fox Times
Sub Product Dimension Processors Time mflops
sgemm() 800 2 × 2 295 4,173
sgemm() 800 4 × 4 121 8,193
sgemm() 1600 2 × 2 1,635 5,010
sgemm() 1600 4 × 4 578 14,173
jki loops 800 2 × 2 980 1,043
jki loops 800 4 × 4 306 3,344
jki loops 1600 2 × 2 7,755 1,056
jki loops 1600 4 × 4 2,103 3,895
92. call mpi_sendrecv_replace(local_B,n_bar*n_bar,mpi_real,&
93. dest, 0,source, 0, &
94. grid%col_comm,status, ierr)
95. end do
96. end subroutine Fox
97. !
98. end program foxmpi
Eight executions were done with foxmpi.f, and the outputs are recorded in
Table 7.6.1. The first four used the sgemm() and the second four used jki loops.
The optimized sgemm() was about four times faster than the jki loops. The
time units were in milliseconds, and the mflops for the larger dimensions were
always the largest.
7.6.5 Exercises
1.
in foxmpi.f.
2. Duplicate the computations for foxmpi.f, and also use a 8 × 8 grid of
processors.
3. Compare the matrix product scheme used in Section 6.5 mmmpi.f with
the Fox algorithm in foxmpi.f.
© 2004 by Chapman & Hall/CRC
Verify for n=6that the matrix product is correct. See lines 23 and 35
Chapter 8
Classical Methods for
Ax = d
The first three sections contain a description of direct methods based on the
Schur complement and domain decomposition. After the first section the co-
e
!cient matrix will be assumed to b e symmetric positive definite (SPD). In
Section 8.3 an MPI code will be studied that illustrates domain decomposition.
Iterative methods based on P-regular splittings and domain decompositions will
be described in the last three sections. Here convergence analysis will be given
via the minimization of the equivalent quadratic functional. An MPI version of
SOR using domain decomposition will be presented in Section 8.5. This chapter
is more analysis-oriented and less application-driven.
8.1 Gauss Elimination
Gauss elimination method, which was introduced in Section 2.2, requires D to
be factored into a product of a lower and upper triangular matrices that have
inverses. This is not always possible, for example, consider the 2 × 2 matrix
where the
d
11
is zero. Then one can interchange the second and first rows
0 1
1 0
¸
0 d
12
d
21
d
22
¸
=
d
21
d
22
0 d
12
¸
=
If d
11
is not zero, then one can use an elementary row operation
1 0
d
21
@d
11
1
¸
d
11
d
12
d
21
d
22
¸
=
d
11
d
12
0 d
22
(d
21
@d
11
)d
12
¸
=
Definition. If there is a permutation matrix S such that S D = OX where O
and X are invertible lower and upper triangular matrices, respectively, then the
matrix
S D is said to have an OX factorization.
313
© 2004 by Chapman & Hall/CRC
314 CHAPTER 8. CLASSICAL METHODS FOR AX = D
Gaussian elimination method for the solution of
D{ = g uses permutations
of the rows and elementary row operations to find the LU factorization so that
S D{ = O(X{) = S g
solve O| = S g and
solve
X{ = |=
Example. Consider the 3 × 3 matrix
5
7
1 2 0
1 2 1
0 1 3
6
8
=
Since the component in the first row and column is not zero, no row interchange
is necessary for the first column. The elementary row operation on the first
column is
5
7
1 0 0
1 1 0
0 0 1
6
8
5
7
1 2 0
1 2 1
0 1 3
6
8
=
5
7
1 2 0
0 0 1
0 1 3
6
8
=
For column two we must interchange rows two and three
5
7
1 0 0
0 0 1
0 1 0
6
8
5
7
1 0 0
1 1 0
0 0 1
6
8
5
7
1 2 0
1 2 1
0 1 3
6
8
=
5
7
1 2 0
0 1 3
0 0 1
6
8
=
Note the first two factors on the left side can be rewritten as
5
7
1 0 0
0 0 1
0 1 0
6
8
5
7
1 0 0
1 1 0
0 0 1
6
8
=
5
7
1 0 0
0 1 0
1 0 1
6
8
5
7
1 0 0
0 0 1
0 1 0
6
8
=
This gives the desired factorization of the matrix
5
7
1 0 0
0 1 0
1 0 1
6
8
5
7
1 0 0
0 0 1
0 1 0
6
8
5
7
1 2 0
1 2 1
0 1 3
6
8
=
5
7
1 2 0
0 1 3
0 0 1
6
8
5
7
1 0 0
0 0 1
0 1 0
6
8
5
7
1 2 0
1 2 1
0 1 3
6
8
=
5
7
1 0 0
0 1 0
1 0 1
6
8
5
7
1 2 0
0 1 3
0 0 1
6
8
=
In order to extend this to q × q matrices, as in Section 2.4 consider just a
2 × 2 block matrix where the diagonal blocks are square but may not have the
same dimension
D =
E H
I F
¸
= (8.1.1)
© 2004 by Chapman & Hall/CRC
8.1. GAUSS ELIMINATION 315
In general
D is q × q with q = n + p, E is n × n, F is p × p, H is n × p and
I is p × n= If E has an inverse, then we can multiply block row one by I E
1
and subtract it from block row two. This is equivalent to multiplication of D
by a block elementary matrix of the form
L
n
0
I E
1
L
p
¸
=
If D{ = g is viewed in block form, then
E H
I F
¸
[
1
[
2
¸
=
G
1
G
2
¸
= (8.1.2)
The above block elementary matrix multiplication gives
E H
0 F I E
1
H
¸
[
1
[
2
¸
=
G
1
G
2
I E
1
G
1
¸
= (8.1.3)
So, if the block upper triangular matrix has an inverse, then this last block
equation can be solved.
The following basic properties of square matrices play an important role in
the solution of (8.1.1). These properties follow directly from the definition of
an inverse matrix.
Theorem 8.1.1 (Basic Matrix Properties) Let B and C be square matrices
that have inverses. Then the following equalities hold:
1.
E
0
0
F
¸
1
=
E
1
0
0
F
1
¸
>
2.
L
n
0
I L
p
¸
1
=
L
n
0
I L
p
¸
>
3.
E
0
I F
¸
=
E
0
0
F
¸
L
n
0
F
1
I L
p
¸
and
4.
E
0
I F
¸
1
=
E
1
0
F
1
I E
1
F
1
¸
=
Definition. Let D have the form in (8.1.1) and E be nonsingular. The Schur
complement of
E in D is F I E
1
H.
Theorem 8.1.2 (Schur Complement Existence) Consider
D as in (8.1.1) and
let
E have an inverse. D has an inverse if and only if the Schur complement of
E in D has an inverse.
© 2004 by Chapman & Hall/CRC
316 CHAPTER 8. CLASSICAL METHODS FOR AX = D
Proof. The proof of the Schur complement theorem is a direct consequence
of using a block elementary row operation to get a zero matrix in the block row
2 and column 1 position
L
n
0
I E
1
L
p
¸
E H
I F
¸
=
E H
0 F I E
1
H
¸
=
Assume that D has an inverse and show the Schur complement must have an
inverse. Since the two matrices on the left side have inverses, the matrix on
the right side has an inverse. Because E has an inverse, the Schur complement
must have an inverse. Conversely,
D may be factored as
E H
I F
¸
=
L
n
0
I E
1
L
p
¸
E H
0 F I E
1
H
¸
=
If both E and the Schur complement have inverses, then both matrices on the
right side have inverses so that D also has an inverse.
The choice of the blocks E and F can play a very important role. Often the
choice of the physical object, which is being modeled, does this. For example
consider the airflow over an aircraft. Here we might partition the aircraft into
wing, rudder, fuselage and "connecting" components. Such partitions of the
physical object or the matrix are called domain decompositions.
Physical problems often should have a unique solution, and so we shall
assume that for
D{ = g, if it has a solution, then the solution is unique. This
means if
D{ = g and Db{ = g> then { = b{= This will be true if the only solution
of D} = 0 is the zero vector } = 0 because for } = { b{
D{ D
b{ = D({ b{) = g g = 0=
Another important consequence of D} = 0 implies } = 0 is that no column of D
can be the zero column vector. This follows from contradiction, if the column
j of
D is a zero vector, then let } be the j unit column vector so that D} = 0
and
} is not a zero vector.
The condition,
D} = 0 implies } = 0> is very helpful in finding the factor-
ization
S D = OX= For example, if D is
2 × 2 and
D} = 0 implies } = 0> then
column one must have at least one nonzero component so that either
S D = X
or D = OX> see the first paragraph of this section. If D has an inverse, then
D} = 0 implies } = 0 and either SD = X or O
1
D = X so that X must have
an inverse. This generalizes via mathematical induction when
D is q × q=
Theorem 8.1.3 (LU Factorization) If D has an inverse so that D} = 0 implies
} = 0> then there exist permutation matrix S> and invertible lower and upper
triangular matrices
O and X, respectively, such that SD = OX=
Proof. We have already proved the q = 2 case. Assume it is true for
any (
q 1) × (q 1) matrix. If D is invertible, then column one must have
some nonzero component, say in row i. So, if the first component is not zero,
© 2004 by Chapman & Hall/CRC
8.1. GAUSS ELIMINATION 317
interchange the first row with the row i. Let
S be the associated permutation
matrix so that
S D =
e h
i F
¸
where
e 6= 0 is 1×1> h is 1×(q1)> i is (q1)×1 and F is (q1)×(q1)= Apply
the Schur complement analysis and the block elementary matrix operation
1 0
ie
1
L
q1
¸
e h
i F
¸
=
e h
0
b
F
¸
where
b
F = F i e
1
h is the (q 1) × (q 1) Schur complement. By the Schur
complement theorem
b
F must have an inverse. Use the inductive assumption to
write
b
F =
b
S
b
O
b
X=
e h
0
b
F
¸
=
e h
0
b
S
b
O
b
X
¸
=
1 0
0
b
S
¸
1 0
0
b
O
¸
e h
0
b
X
¸
=
Since
b
S is a permutation matrix,
1 0
0
b
S
¸
1 0
ie
1
L
q1
¸
S D =
1 0
0
b
O
¸
e h
0
b
X
¸
=
Note
1 0
0
b
S
¸
1 0
ie
1
L
q1
¸
=
1 0
b
S (ie
1
)
b
S
¸
=
1 0
b
S (ie
1
) L
q1
¸
1 0
0
b
S
¸
=
Then
1 0
b
S (ie
1
) L
q1
¸
1 0
0
b
S
¸
S D =
1 0
0
b
O
¸
e h
0
b
X
¸
=
Finally, multiply by the inverse of the left factor on the left side to get the
desired factorization
1 0
0
b
S
¸
S D =
1 0
b
S ie
1
L
q1
¸
1 0
0
b
O
¸
e h
0
b
X
¸
=
1 0
b
S ie
1
b
O
¸
e h
0
b
X
¸
=
In order to avoid significant roundo errors due to small diagonal compo-
nents, the row interchanges can be done by choosing the row with the largest
© 2004 by Chapman & Hall/CRC
318 CHAPTER 8. CLASSICAL METHODS FOR AX = D
possible nonzero component. If the matrix is symmetric positive definite, then
the diagonals will be p ositive and the row interchanges may not be necessary.
In either case one should give careful consideration to u sing the subroutines in
LAPACK [1].
8.1.1 Exercises
1. Consider the 3 × 3 example. Use the Schur complement as in (8.1.3) to
solve
D{ = g = [1 2 3]
W
where E = [1]=
2. Consider the 3 × 3 example. Identify the steps in the existence theorem
for the factorization of the matrix.
3. In the proof of the existence theorem for the factorization of the matrix,
write the factorization in component form using the Schur complement
b
F = F (i@e)h = [f
lm
(i
l
@e)h
m
]
where
h
m
is mth component of the 1 × (q 1) array h= Write this as either the ij
or ji version and explain why these could be described as the row and column
versions, respectively.
4. Assume
E is an invertible n × n matrix.
(a). Verify the following for
b
F = F I E
1
I
L
n
0
I E
1
L
p
¸
E H
I F
¸
L
n
E
1
H
0 L
p
¸
=
E
0
0
b
F
¸
=
(b). Use this to show D has an i nverse if and only if
b
F has an inverse.
5. Assume
D is an q × q matrix and prove the following are equivalent:
(i).
D has an inverse,
(ii). there exist permu tation matrix
S> and
invertible lower and upper triangular matrices
O and X,
respectively, such that
S D = OX and
(iii).
D} = 0 implies } = 0=
8.2 Symmetric Positive Definite Matrices
In this section we will restrict the matrices to symmetric positive definite ma-
trices. Although this restriction may seem a little severe, there are a num-
ber of important applications, which include some classes of partial di
erential
equations and some classes of least squares problems. The advantage of this
restriction is that the number of operations to do Gaussian elimination can be
cut in half.
Definition. Let
D be an q × q real matrix. D is a real symmetric positive
definite matrix (SPD) if and only if
D = D
W
and for all { 6= 0> {
W
D{ A 0.
Examples.
© 2004 by Chapman & Hall/CRC
8.2. SYMMETRIC POSITIVE DEFINITE MATRICES 319
1. Consider the 2 × 2 matrix
2 1
1 2
¸
and note
{
W
D{ = {
2
1
+ ({
1
{
2
)
2
+ {
2
2
A 0=
A similar q × q matrix is positive definite
5
9
9
9
9
7
2 1
1 2
.
.
.
.
.
.
.
.
.
1
1 2
6
:
:
:
:
8
=
2. The matrix D for which D = D
W
with d
ll
A
P
m6=l
|d
lm
| is positive definite.
The symmetry implies that the matrix is also column strictly diagonally domi-
nant, that is, d
mm
A
P
l6=m
|d
lm
˙
|. Now use this and the inequality |
de|
1
2
(d
2
+e
2
)
to show for all
{ 6= 0> {
W
D{ A 0.
3. Consider the normal equation from the least squares problem where
D is
p × q where p A q. Assume A has full column rank (D{ = 0 implies { = 0),
then the normal equation
D
W
D{ = D
W
g is equivalent to finding the least squares
solution of
D{ = g. Here D
W
D is SPD because if { 6= 0 , then D{ 6= 0, and
{
W
(D
W
D){ = (D{)
W
(D{) A 0=
Theorem 8.2.1 (Basic Properties of SPD Matrices) If A is an q × q SPD
matrix, then
1. The diagonal components of
D are positive, d
ll
A 0,
2. If
D =
E I
W
I F
¸
> then E and F are SPD,
3.
D{ = 0 implies { = 0 so that D has an inverse and
4. If
V is p×q with p q and has full column rank, then V
W
DV is positive
definite.
Proof. 1. Choose
{ = h
l
, the unit vector with 1 in component i so that
{
W
D{ = d
ll
A 0=
2. Choose { =
[
1
0
¸
so that
{
W
D{ = [
W
1
E[
1
A 0=
3. Let { 6= 0 so that by the positive definite assumption {
W
D{ A 0 and
D{ 6= 0. This is equivalent to item 3 and the existence of an inverse matrix.
4. Let
{ 6= 0 so that by the full rank assumption on V V{ 6= 0= By the
positive definite assumption on
D
(V{)
W
D(V{) = {
W
(V
W
DV){ A 0=
© 2004 by Chapman & Hall/CRC
320 CHAPTER 8. CLASSICAL METHODS FOR AX = D
The next theorem uses the Schur complement to give a characterization of
the block 2 × 2 SPD matrix
D =
E I
W
I F
¸
= (8.2.1)
Theorem 8.2.2 (Schur Complement Characterization) Let
D as in (8.2.1) be
symmetric.
D is SPD if and only if E and the Schur complement of E in D>
b
F = F I E
1
I
W
> are SPD.
Proof. Assume D is SPD so that E is also SPD and has an inverse. Then
one can use block row and column elementary operations to show
L
n
0
I E
1
L
p
¸
E I
W
I F
¸
L
n
E
1
I
W
0 L
p
¸
=
E
0
0
b
F
¸
=
Since E is SPD,
¡
E
1
¢
W
=
¡
E
W
¢
1
= E
1
and thus
V =
L
n
E
1
I
W
0 L
p
¸
and
V
W
=
L
n
E
1
I
W
0 L
p
¸
W
=
L
n
0
I E
1
L
p
¸
=
Then
V
W
DV =
E
0
0
b
F
¸
=
Since V has an inverse, it has full rank and E and
b
F must be SPD. The converse
is also true by reversing the above argument.
Example. Consider the 3 × 3 matrix
D =
5
7
2 1 0
1 2 1
0
1 2
6
8
=
The first elementary row operation is
5
7
1 0 0
1
@2 1 0
0 0 1
6
8
5
7
2 1 0
1 2 1
0
1 2
6
8
=
5
7
2 1 0
0 3
@2 1
0
1 2
6
8
=
The Schur complement of the 1 × 1 matrix E = [2] is
b
F = F I E
1
I
W
=
3@2 1
1 2
¸
=
It clearly is SPD and so one can do another elementary row operation, but now
on the second column
5
7
1 0 0
0 1 0
0 2@3 1
6
8
5
7
2 1 0
0 3
@2 1
0
1 2
6
8
=
5
7
2 1 0
0 3
@2 1
0 0 4
@3
6
8
=
© 2004 by Chapman & Hall/CRC