White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI

Подождите немного. Документ загружается.

300 CHAPTER 7. MESSAGE PASSING INTERFACE

12. call mpi_init(ierr)

13. call mpi_comm_rank(mpi_comm_world,my_rank,ierr)

14. call mpi_comm_size(mpi_comm_world,p,ierr)

15.! Pro cessor 0 packs and broadcasts the four number.

16. if (my_rank.eq.0) then

17. a = exp(1.)

18. b = 4*atan(1.)

19. c = 1

20. d = 186000

21. location = 0

22. call mpi_pack(a,1,mpi_real,numbers,100,location,&

mpi_comm_world, ierr)

23. call mpi_pack(b,1,mpi_real,numbers,100,location,&

mpi_comm_world, ierr)

24. call mpi_pack(c,1,mpi_integer,numbers,100,location,&

mpi_comm_world, ierr)

25. call mpi_pack(d,1,mpi_integer,numbers,100,location,&

mpi_comm_world, ierr)

26. call mpi_bcast(numbers,100,mpi_packed,0,&

mpi_comm_world,ierr)

27. else

28. call mpi_bcast(numbers,100,mpi_packed,0,&

mpi_comm_world,ierr)

29.! Each processor unpacks the numbers.

30. location = 0

31. call mpi_unpack(numbers,100,location,a,1,mpi_real,&

mpi_comm_world, ierr)

32. call mpi_unpack(numbers,100,location,b,1,mpi_real,&

mpi_comm_world, ierr)

33. call mpi_unpack(numbers,100,location,c,1,mpi_integer,&

mpi_comm_world, ierr)

34. call mpi_unpack(numbers,100,location,d,1,mpi_integer,&

mpi_comm_world, ierr)

35. end if

36.! Each processor prints the numbers.

37. print*,’my_rank =’,my_rank, ’a = ’,a

38. print*,’my_rank =’,my_rank, ’b = ’,b

39. print*,’my_rank =’,my_rank, ’c = ’,c

40. print*,’my_rank =’,my_rank, ’d = ’,d

41. call mpi_ﬁnalize(ierr)

42. end program packmpi

my_rank = 0 a = 2.718281746

my_rank = 0 b = 3.141592741

7.5. COMMUNICATORS 301

my_rank = 0 c = 1

my_rank = 0 d = 186000

my_rank = 1 a = 2.718281746

my_rank = 1 b = 3.141592741

my_rank = 1 c = 1

my_rank = 1 d = 186000

my_rank = 2 a = 2.718281746

my_rank = 2 b = 3.141592741

my_rank = 2 c = 1

my_rank = 2 d = 186000

my_rank = 3 a = 2.718281746

my_rank = 3 b = 3.141592741

my_rank = 3 c = 1

my_rank = 3 d = 186000

7.4.5 Exercises

1. Duplicate the calculations for countmpi.f and experiment with dierent

size arrays and numbers of processors.

2. Duplicate the calculations for dertypempi.f and experiment with di

erent

data types.

3. Duplicate the calculations for packmpi.f and experiment with di

erent

numbers of processors and di

erent size vectors.

4. Consider a one dimensional array that has many nonzero numbers. Use

mpi_pack() and mpi_unpack() to communicate the nonzero entries in the array.

5. Repeat exercise 4 for a two dimensional array.

7.5 Communicators

7.5.1 Intro duction

The generic communicator that has been used in all of the previous codes is

called mpi_comm_world. It is the set of all p processors, and all processors

can communicate with each other. When collective subroutines are called and

the communicator is mpi_comm_world, then the data is to be communicated

among the other p-1 processors. In many applications it may not be necessary

to communicate with all other processors. Two similar examples were given

in the Section 6.6, where the two space dimension heat and pollutant models

are considered. In these cases each processor is associated with a horizontal

portion of space, and each processor is required to exchange information with

the adjacent processors.

302 CHAPTER 7. MESSAGE PASSING INTERFACE

7.5.2 A Grid Communicator

In this and the next section grid communicators will be used to do matrix-

matrix products. In order to motivate this discussion, consider a block 3 × 3

matrix times a block 3 × 1 matrix where the blocks are q × q

[

+ D

[

+ D

[

+ D

[

+ D

[

+ D

[

+ D

[

Consider s = 9 processors and associate then with a 3 × 3 grid. Assume the

matrices

are stored on grid processor lm= Then the 9 matrix pro ducts D

[

could be done concurrently. The overall process is as follows:

broadcast

[

to column m of processors,

in parallel compute the matrix products

[

and

sum the products in row

l of processors.

We wish to use collective subroutines restricted to either columns or rows of this

grid of processors. In particular, start with the generic communicator, and then

create a two dimension grid communicator as well as three row and three column

subgrid communicators. The MPI subroutines mpi_cart_create(), mpi_cart_coords()

and mpi_cart_sub() will help us do this. The subroutine Setup_grid(grid) uses

these three MPI subroutines, and it is used in gridcommpi.f, and in foxmpi.f of

the next section.

The subroutine mpi_cart_create() will generate a

g = glp dimensional

grid communicator from

s = t

processors in the original communicator called

mpi_comm. The

glpvl}h = t with periodic parameter in each dimension is

set equal to TRUE, and the numbering of the pro cessor in each row or column

of processor will begin at 0 and end at

t  1. The new communicator is called

grid_comm.

call mpi_cart_create(mpi_comm_world, dim,&

dimsize, periods, .TRUE. , grid_comm, ierr)

mpi_comm integer

dim integer

dimsize integer(*)

p eriods logical(*)

reorder logical

grid_comm integer

ierr logical

The subroutine mpi_cart_coords() asso ciates with each processor, given by

grid_my_rank, in grid_comm, a grid_row = coordinates(0) or grid_col =

coordinates(1) for

glp = 2.

call mpi_cart_coords(grid_comm, grid_my_rank, 2,&

co ordinates, ierr )

7.5. COMMUNICATORS 303

grid_comm integer

grid_my_rank integer

dim integer

coordinates integer(*)

ierr logical

Subcommunicators can easily be formed by a call to mpi_cart_sub(). The

subcommunicators are associated with the grid row or grid columns of proces-

sors for dim = 2.

call mpi_cart_sub(grid_comm, vary_coords, &

sub_comm, ierr)

grid_comm integer

vary_coords logical(*)

sub_comm integer

ierr logical

7.5.3 Illustration gridcommpi.f

First, we examine the subroutine Setup_grid(grid) in lines 56-91. This sub-

routine and subsequent subroutines on communicators are Fortran 9x varia-

tions of those in Pacheco [21]. This is a two dimensional grid and we have

assumed

s = t

. The parameter grid is of type GRID_INFO_TYPE as de-

ﬁned in lines 5-14. The integer array dimension and logical array periods are

deﬁned in lines 70-73. The call to mpi_cart_create() is done in line 74 where

a grid%comm is deﬁned. In line 76 grid%my_rank is identiﬁed for the commu-

nicator grid%comm. Lines 79-80 identify the grid%my_row and grid%my_col.

In lines 84 and 88 mpi_cart_sub() deﬁne the communicators grid%row_comm

and grid%col_comm.

Second, the main part of gridcommpi.f simply deﬁnes a 6 × 6 array and uses

s = 3

processors so that the array can be deﬁned by nine processors as given

in lines 26-32. The local arrays

D are 2 × 2, and there is a version on each of

the nine processors. After line 32 the 6 × 6 array, which is distributed over the

grid communicator, is

1 1 2 2 3 3

2 2 3 3 4 4

3 3 4 4 5 5

In line 48 mpi_bcast() from column processors 1 (corresponds to the second

block column in the above matrix) to the other processors in row_comm. This

304 CHAPTER 7. MESSAGE PASSING INTERFACE

means the new distribution of the matrix will be

2 2 2 2 2 2

3 3 3 3 3 3

4 4 4 4 4 4

The output from the print command in lines 50-53 veriﬁes this.

MPI/Fortran 9x gridcommpi.f

1. program gridcommpi

2.! Illustrates grid commu nicators.

3. include ’mpif.h’

4. IMPLICIT NONE

5. type GRID_INFO_TYPE

6. integer p ! total number of processes.

7. integer comm ! communicator for the entire grid.

8. integer row_comm ! communicator for my row.

9. integer col_comm ! communicator for my col.

10. integer q ! order of grid.

11. integer my_row ! my row number.

12. integer my_col ! my column number.

13. integer my_rank ! my rank in the grid communicator.

14. end type GRID_INFO_TYPE

15. TYPE (GRID_INFO_TYPE) :: grid_info

16. integer :: my_rank, ierr

17. integer, allocatable, dimension(:,:) :: A,B,C

18. integer :: i,j,k,n, n_bar

19. call mpi_init(ierr)

20. call Setup_grid(grid_info)

21. call mpi_comm_rank(mpi_comm_world, my_rank, ierr)

22. if (my_rank == 0) then

23. n=6

24. endif

25. call mpi_bcast(n,1,mpi_integer, 0, mpi_comm_world, ierr)

26. n_bar = n/(grid_info%q)

27.! Allocate local storage for local matrix.

28. allocate( A(n_bar,n_bar) )

29. allocate( B(n_bar,n_bar) )

30. allocate( C(n_bar,n_bar) )

31. A = 1 + grid_info%my_row + grid_info%my_col

32. B = 1 - grid_info%my_row - grid_info%my_col

33. if (my_rank == 0) then

34. print*,’n = ’,n,’n_bar = ’,n_bar,&

7.5. COMMUNICATORS 305

35. ’grid%p = ’,grid_info%p, ’grid%q = ’,grid_info%q

36. end if

37. print*, ’my_rank = ’,my_rank,&

38. ’grid_info%my_row = ’,grid_info%my_row,&

39. ’grid_info%my_col = ’,grid_info%my_col

40. call mpi_barrier(mpi_comm_world, ierr)

41. print*, ’grid_info%my_row =’,grid_info%my_row,&

42. ’grid_info%my_col =’,grid_info%my_col,&

43. ’A = ’,A(1,:),&

44. ’ ; ’,A(2,:)

45.! Uses mpi_bcast to send and receive parts of the array, A,

46.! to the processors in grid_info%row_com, which was deﬁned

47.! in the call to the subroutine Setup_grid(grid_info).

48. call mpi_bcast(A,n_bar*n_bar,mpi_integer,&

49. 1, grid_info%row_comm, ierr)

50. print*, ’grid_info%my_row =’,grid_info%my_row,&

51. ’grid_info%my_col =’,grid_info%my_col,&

52. ’ new_A = ’,A(1,:),&

53. ’ ; ’,A(2,:)

54. call mpi_ﬁnalize(ierr)

55. contains

56. subroutine Setup_grid(grid)

57.! This subroutine deﬁnes a 2D grid communicator.

58.! And for each grid row and grid column additional

59.! communicators are deﬁned.

60. TYPE (GRID_INFO_TYPE), intent(inout) :: grid

61. integer old_rank

62. integer dimensions(0:1)

63. logical periods(0:1)

64. integer coordinates(0:1)

65. logical varying_coords(0:1)

66. integer ierr

67. call mpi_comm_size(mpi_comm_world, grid%p, ierr)

68. call mpi_comm_rank(mpi_comm_world, old_rank, ierr )

69. grid%q = int(sqrt(dble(grid%p)))

70. dimensions(0) = grid%q

71. dimensions(1) = grid%q

72. periods(0) = .TRUE.

73. periods(1) = .TRUE.

74. call mpi_cart_create(mpi_comm_world, 2,&

75. dimensions, periods, .TRUE. , grid%comm, ierr)

76. call mpi_comm_rank (grid%comm, grid%my_rank, ierr )

77. call mpi_cart_coords(grid%comm, grid%my_rank, 2,&

78. coordinates, ierr )

306 CHAPTER 7. MESSAGE PASSING INTERFACE

79. grid%my_row = coordinates(0)

80. grid%my_col = coordinates(1)

81.! Set up row and column communicators.

82. varying_coords(0) = .FALSE.

83. varying_coords(1) = .TRUE.

84. call mpi_cart_sub(grid%comm,varying_coords,&

85. grid%row_comm,ierr)

86. varying_coords(0) = .TRUE.

87. varying_coords(1) = .FALSE.

88. call mpi_cart_sub(grid%comm,varying_coords,&

89. grid%col_comm,ierr)

90. end subroutine Setup_grid

91. end program gridcommpi

n = 6 n_bar = 2 grid%p = 9 grid%q = 3

my_rank = 0 grid_info%my_row = 2 grid_info%my_col = 2

my_rank = 1 grid_info%my_row = 2 grid_info%my_col = 1

my_rank = 2 grid_info%my_row = 2 grid_info%my_col = 0

my_rank = 3 grid_info%my_row = 1 grid_info%my_col = 0

my_rank = 4 grid_info%my_row = 1 grid_info%my_col = 2

my_rank = 5 grid_info%my_row = 0 grid_info%my_col = 2

my_rank = 6 grid_info%my_row = 0 grid_info%my_col = 1

my_rank = 7 grid_info%my_row = 0 grid_info%my_col = 0

my_rank = 8 grid_info%my_row = 1 grid_info%my_col = 1

grid_info%my_row = 0 grid_info%my_col = 0 A = 1 1 ; 1 1

grid_info%my_row = 1 grid_info%my_col = 0 A = 2 2 ; 2 2

grid_info%my_row = 2 grid_info%my_col = 0 A = 3 3 ; 3 3

grid_info%my_row = 0 grid_info%my_col = 1 A = 2 2 ; 2 2

grid_info%my_row = 1 grid_info%my_col = 1 A = 3 3 ; 3 3

grid_info%my_row = 2 grid_info%my_col = 1 A = 4 4 ; 4 4

grid_info%my_row = 0 grid_info%my_col = 2 A = 3 3 ; 3 3

grid_info%my_row = 1 grid_info%my_col = 2 A = 4 4 ; 4 4

grid_info%my_row = 2 grid_info%my_col = 2 A = 5 5 ; 5 5

grid_info%my_row = 0 grid_info%my_col = 0 new_A = 2 2 ; 2 2

grid_info%my_row = 1 grid_info%my_col = 0 new_A = 3 3 ; 3 3

grid_info%my_row = 2 grid_info%my_col = 0 new_A = 4 4 ; 4 4

grid_info%my_row = 0 grid_info%my_col = 1 new_A = 2 2 ; 2 2

grid_info%my_row = 1 grid_info%my_col = 1 new_A = 3 3 ; 3 3

grid_info%my_row = 2 grid_info%my_col = 1 new_A = 4 4 ; 4 4

grid_info%my_row = 0 grid_info%my_col = 2 new_A = 2 2 ; 2 2

grid_info%my_row = 1 grid_info%my_col = 2 new_A = 3 3 ; 3 3

grid_info%my_row = 2 grid_info%my_col = 2 new_A = 4 4 ; 4 4

7.6. FOX ALGORITHM FOR AB 307

7.5.4 Exercises

1. Duplicate the computations for gridcommpi.f. Change mpi_bcast() to

mpi_bcast(A,n_bar*n_bar,mpi_real,x,grid_info%row_comm,ierr) where x is

0 and 2. Explain the outputs.

2. In gridcommpi.f change the communicator from row_comm to col_comm.by

using mpi_bcast(A,n_bar*n_bar,mpi_real,x,grid_info%col_comm,ierr) where

x is 0, 1 and 2. Explain the outputs.

7.6 Fox Algorithm for AB

7.6.1 Intro duction

In this section the block matrix-matrix product DE will be done where D and

E are both t × t block matrices. The number of processors used to do this

will be s = t

, and the grid communicator that was deﬁned in the subroutine

Setup_grid() will be used. Fox’s algorithm follows a similar pattern as in

in the previous section where D is 3 × 3 and [ is 3 × 1= The numbering of the

block rows and columns of the matrices will start at 0 and end at

t  1.

7.6.2 Matrix-Matrix Product

The classical deﬁnition of matrix product F = DE is block lm of F equals block

row

l of D times block column m of E

t1

n=0

The summation can be done in any order, and the matrix products for a ﬁxed

lm can be done by the grid processor lm= If the matrices D

, E

, and F

are stored on the grid processor lm, then the challenge is to communicate the

required matrices to the grid processor

lm=

7.6.3 Parallel Fox Algorithm

In order to motivate the Fox algorithm, consider the block 3 × 3 case

Assume that pro cessor lm in the grid communicator has stored the matrices D

and E

. Consider the computation of the second block row of F, which can be

308 CHAPTER 7. MESSAGE PASSING INTERFACE

reordered as follows

= D

+ D

= D

+ D

= D

+ D

Grid processor 1m can compute the ﬁrst term on the right side if the matrix

has been broadcast to grid processor 1m. In order for grid processor 1m to

compute the second matrix product on the right side, the matrix E

must be

replaced by

, and matrix D

must be broadcast to the grid processor 1m.

The last step is for t  1 = 2 where the matrix E

must be replaced by E

and the matrix

must be broadcast to the grid processor 1m= For the t × t

block matrices there are t matrix products for each grid processor, and this can

be done in a loop whose index is step as in the following algorithm.

Fox Algorithm for the Matrix Product

F = F + DE

= s

1@2

, vrxufh = (prg(l + 1> t)> m), ghvw = (prg(l  1> t)> m)

concurrently with 0  l> m  t  1

for

vwhs = 0> t  1

n_edu = prg(l + vwhs> t)

broadcast

l>n_edu

to grid row l of processors

= F

+ D

l>n_edu

n_edu>m

send E

n_edu>m

to processor dest

receive

nn>m

from source where nn = prg(n_edu + 1> t)

endloop.

7.6.4 Illustration foxmpi.f

This implementation of the Fox algorithm is a variation on that given by

Pacheco [21], but here Fortran 9x is used, and the matrix products for the

submatrices are done either by a call to the BLAS3 subroutine sgemm() or

by the jki loops. The input to this is given in lines 1-33, the call to the fox

subroutine is in line 34, and the output in given in lines 36-42. The subroutine

Setup_grid(grid) is the same as listed in the previous section. The Fox sub-

routine is listed in lines 48-96. The step loop of the Fox algorithm is executed

in lines 60-95. The matrix pro ducts may be done by either sgemm() or the jki

loops, which are listed here as commented out of execution. The broadcast of

the matrix

l>n_edu

over the grid_row communicator is done in lines 63-64 and

79-80; note how one stores local_A from bcast_root in temp_A so as not to

overwrite local_A in destination. Then the matrix products

l>n_edu

n_edu>m

are done in lines 66-67 and 81-82. The subroutine mpi_sendrecv_replace() in

lines 92-94 is used to communicate E

n_edu>m

within the grid_col communicator.

MPI/Fortran 9x Code foxmpi.f

1. program foxmpi

7.6. FOX ALGORITHM FOR AB 309

2. include ’mpif.h’

3. IMPLICIT NONE

4. type GRID_INFO_TYPE

5. integer p ! total number of processes.

6. integer comm ! communicator for the entire grid.

7. integer row_comm ! communicator for my row.

8. integer col_comm ! communicator for my col.

9. integer q ! order of grid.

10. integer my_row ! my row number.

11. integer my_col ! my column number.

12. integer my_rank ! my rank in the grid communicator.

13. end type GRID_INFO_TYPE

14. TYPE (GRID_INFO_TYPE) :: grid_info

15. integer :: my_rank, ierr

16. real, allocatable, dimension(:,:) :: A,B,C

17. integer :: i,j,k,n, n_bar

18. real:: mﬂops,t1,t2,timef

19. call mpi_init(ierr)

20. call Setup_grid(grid_info)

21. call mpi_comm_rank(mpi_comm_world, my_rank, ierr)

22. if (my_rank == 0) then

23. n = 800 !n = 6

24. t1 = timef()

25. endif

26. call mpi_bcast(n,1,mpi_integer, 0, mpi_comm_world,ierr)

27. n_bar = n/(grid_info%q)

28. ! Allocate storage for local matrix.

29. allocate( A(n_bar,n_bar) )

30. allocate( B(n_bar,n_bar) )

31. allocate( C(n_bar,n_bar) )

32. A = 1.0 + grid_info%my_row + grid_info%my_col

33. B = 1.0 - grid_info%my_row - grid_info%my_col

34. call Fox(n,grid_info,A,B,C,n_bar)

35. ! print*,grid_info%my_row, grid_info%my_col, ’C = ’,C

36. if (my_rank == 0) then

37. t2 = timef()

38. print*,t2

39. print*,n,n_bar,grid_info%q

40. mﬂops = (2*n*n*n)*.001/t2

41. print*, mﬂops

42. endif

43. call mpi_ﬁnalize(ierr)

44. contains

45. !

46. !subroutine Setup_grid . . . .see chapter 7.5 and gridcommpi.f