Hands on Intro to parallel computing in FORTRAN

MPI FORTRAN

1. Hands on intro to
parallel computing in
Fortran
time < 3 minutes

2. Problem
Simple algebraic system
•One dimensional steady state heat transfer
or
•Discretize on 1 dimensional grid
•At point i,
∇⋅k(∇T)=0k
∂
2
T
∂x
2
=0
T
i+1−2T
i+T
i−1
δx
2
=0
T
h
T
c
i=0 ii-1 i+1
x
δxδx

3. Single domain
System of Equations
•For i=1 to N

•The equivalent system
T
i−1
−2T
i
+T
i+1
δx
2
−2100...0
1−210...0
0⋱⋱⋱...0
0...01−21
0...001−2
N×N
T
1
T
2
⋮
T
n−1
T
n
N×1
=
−T
h
0
⋮
0
−T
c
i=0 ii-1 i+1
δxδx
T
h T
c
i=N+1

4. Single domain
Solution of the System
•Solve directly!
•Or use Jacobi or Gauss Seidel method
•For a given system of linear equations

•The iterative solution is
a
11
x
1
+a
12
x
2
+a
13
x
3
=b
1
a
21
x
1
+a
22
x
2
+a
23
x
3
=b
2
a
31
x
1
+a
32
x
2
+a
33
x
3
=b
3
x
new
1
=x
old
1
+
1
a
11
[
b
1−(a
11x
old
1
+a
12x
old
2
+a
13x
old
3)]
x
new
2
=x
old
2
+
1
a
21
[
b
2
−(a
21
x
old
1
+a
22
x
old
2
+a
23
x
old
3)]
x
new
3
=x
old
3
+
1
a
31
[
b
3
−(a
31
x
old
1
+a
32
x
old
2
+a
33
x
old
3)]
a
11a
12a
13
a
21a
22a
23
a
31
a
32
a
33
x
1
x
2
x
3
=
b
1
b
2
b
3
x
1
x
2
x
3

5. <
Parallel Computation
MPI
i=0 ii-1 i+1
δx
T
h T
c
i=N+1
ii-1 i+1
i=0 i=n+1
Processor 1 Processor 2 Processor M
ii-1 i+1
i=0 i=n+1
n=
N
M
procs
MPI Send
MPI Recv
n+2 nodes

Aqeel Ahmed

6. Parallel Computation
Fortran and MPI
1 program laplacian_1d_parallel
2 use mpi_f08 ! call the MPI library
3 implicit none ! always declare variables (no auto type deduction)
4
5 ! declare variables
6 integer :: points_total, npoints, i, j, k, l, iter_max
7 real, dimension(:), allocatable :: T, T_full
8 real :: T_initial, Tr, Tl, omega, size, Temp, R
9 integer myrank, nranks, tag, ierror
10 type(MPI_Status) :: istatus
11
12 ! initialize MPI
13 call MPI_INIT(ierror)
14 call MPI_COMM_SIZE(MPI_COMM_WORLD, nranks, ierror) ! get no. of procs
15 call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierror) ! get rank
16
17 ! total number of points
18 points_total = 2400
19 ! number of points per processor
20 npoints = int(points_total/nranks)
21 ! allocate T with 2 ghost points T(1) and T(npoints+2)
22 allocate(T(npoints+2))
23 if (myrank == 0) then
24 allocate(T_full(npoints*nranks)) ! this is for gather only
25 end if
26
27 ! initial conditions and parameters
28 Tl = 25;Tr = 18;T_initial = 20;omega = 1.5;Temp = 0.0;iter_max = 500000
29 T = T_initial
30 R = 0.0
31
32 ! set boundary conditions
33 if (myrank == 0) then
34 T(1) = Tl
35 end if
36 if (myrank == nranks-1) then
37 T(npoints+2) = Tr
38 end if
39
MPI
specific
variables
Variables
declarations
Domain
discretisation
parameters
Problem
inputs
Boundary
conditions

7. Parallel Computation
Fortran and MPI
40 ! iterative solution of the system
41 do k = 1, iter_max
42 ! trasfer data at processors boundaries
43 do i = 0,nranks-2
44 if (myrank==i) then ! send the 2nd last point to next processor 1st point
45 Temp = T(npoints+1)
46 call MPI_SEND(Temp, 1, MPI_REAL, myrank+1, i, MPI_COMM_WORLD, ierror)
47 elseif (myrank==i+1) then ! receive the 2nd last point from the previous proc as 1st point
48 call MPI_RECV(Temp, 1, MPI_REAL, myrank-1, i, MPI_COMM_WORLD, istatus, ierror)
49 T(1) = Temp
50 end if
51 end do
52
53 do l = 1,nranks-1
54 ! send the second point T(2) of the current proc to previous proc last point T(npoints+2)
55 if (myrank==l) then ! send the second point to previous procs
56 Temp = T(2)
57 call MPI_SEND(Temp, 1, MPI_REAL, myrank-1, l, MPI_COMM_WORLD, ierror)
58 elseif (myrank==l-1) then ! receive the second point from the next procs as last point
59 call MPI_RECV(Temp, 1, MPI_REAL, myrank+1, l, MPI_COMM_WORLD, istatus, ierror)
60 T(npoints+2) = Temp
61 end if
62 end do
63
64 ! solution using Gauss-Seidel method with Successive Over Relaxation (SOR)
65 do i = 2,npoints+1
66 R = 0.5*omega*(T(i+1) - 2*T(i) + T(i-1))
67 T(i) = T(i) + R
68 end do
69
70 ! set boundary conditions
71 if (myrank == 0) then
72 T(1) = Tl
73 end if
74 if (myrank == nranks-1) then
75 T(npoints+2) = Tr
76 end if
77 end do
78
79 ! put T from all procs to root (master) in T_full
80 call MPI_Gather(T(2:npoints+1), npoints, MPI_REAL, T_full, npoints, &
81 MPI_REAL, 0, MPI_COMM_WORLD, ierror)
82
83 call MPI_FINALIZE(ierror)
MPI
communications
Iterative solution