Commit bf302f59 authored by Ciarán Ó Rourke's avatar Ciarán Ó Rourke
Browse files

Merge branch '5-mpi-non-blocking-and-collective-fortran-exercises' into 'master'

Resolve "MPI: non-blocking and collective fortran exercises"

Closes #5

See merge request training/sohpc-training-2021!5
parents aaf12ff8 41f0a37e
......@@ -2,6 +2,6 @@
## Configuration
Parameters can be edited in the [param.h](./params.h) file. `nx`, `ny`, and `nt`
are the size of the grid in the x and y dimensions and the number of timesteps
respectively.
Parameters can be edited in the [param.h](./params.h) file or the main Fortran
file for C and Fortran respectively. `nx`, `ny`, and `nt` are the size of the
grid in the x and y dimensions and the number of timesteps respectively.
! initialise with randomly with even distribution
subroutine randomise(u)
use params
implicit none
real (kind=8), intent(out) :: u(grid_size)
real :: rvec(grid_size)
call random_number(rvec)
u = rvec*10.0
return
end subroutine randomise
subroutine print_local_grid(u, nx, ny)
use params
implicit none
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
integer (kind=4) :: i, j
! don't print halo regions
do i = 1,ny
do j = 0, nx-1
write(6,'(f8.3,1x)',advance='NO') u(i * nx + j+1)
end do
write(6,*)
end do
return
end subroutine print_local_grid
subroutine print_grid(u, nx, ny, rank, unisize)
use mpi
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
integer :: ierror
integer (kind=4) :: i
do i = 0,unisize - 1
if (rank .eq. i) then
call print_local_grid(u, nx, ny)
endif
call MPI_Barrier(MPI_COMM_WORLD, ierror)
end do
! last process prints an extra line
if (rank .eq. i) then
call print_local_grid(u, nx, ny)
write(6,*)
endif
call MPI_Barrier(MPI_COMM_WORLD,ierror)
return
end subroutine print_grid
subroutine write_local_grid_to_file(u, nx, ny)
use params
implicit none
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
integer (kind=4) :: i, j
do i = 1,ny
do j = 0, nx-1
write(unitnum,'(f8.3,1x)',advance='NO') u(i * nx + j+1)
end do
write(unitnum,*)
end do
return
end subroutine write_local_grid_to_file
subroutine write_grid_to_file(filename, u, nx, ny, rank, unisize)
use mpi
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
character (len=14), intent(in) :: filename
integer :: ierror
integer (kind=4) :: i
! first process truncates file to size 0 on open
if (rank .eq. 0) then
open(unit=unitnum,file=filename,access='sequential',form='FORMATTED',status='UNKNOWN',&
action='WRITE',iostat=ierror)
if (ierror .eq. 0) then
call write_local_grid_to_file(u, nx, ny)
close(unit=unitnum,status='KEEP')
endif
endif
call MPI_Barrier(MPI_COMM_WORLD, ierror)
do i = 1,unisize-1
if (rank .eq. i) then
open(unit=unitnum,file=filename,access='sequential',form='FORMATTED',status='OLD',&
action='WRITE',position='APPEND',iostat=ierror)
if (ierror .eq. 0) then
call write_local_grid_to_file(u, nx, ny)
close(unit=unitnum,status='KEEP')
endif
endif
call MPI_Barrier(MPI_COMM_WORLD, ierror)
end do
return
end subroutine write_grid_to_file
module params
integer (kind=4), parameter :: nxg = 10, nyg = 13, unitnum = 2
integer (kind=4), parameter :: nt = 1000 ! number of time-steps to evolve the system
integer (kind=4) :: grid_size
real (kind=8), parameter :: a = 1.0 ! diffusion constant
real (kind=8) :: dx, dy, dx2, dy2, dt
contains
subroutine initialise()
dx = 1.0 / real(nxg)
dy = 1.0 / real(nyg)
dx2 = dx * dx
dy2 = dy * dy
dt = dx2*dy2 / (2.0*a*(dx2 + dy2))
end subroutine initialise
end module params
module funcs
contains
subroutine halo_exchange(u, nx, ny, rank, unisize)
use mpi
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(inout) :: u(grid_size)
integer :: ierror, up_neighbour, down_neighbour
real (kind=8) :: buf(nx)
up_neighbour = modulo((rank - 1 + unisize), unisize)
down_neighbour = modulo((rank + 1), unisize)
!
! checkerboard scheme to avoid deadlocks
! note: no longer required for non-blocking
!
if (modulo(rank,2) .eq. 0) then
! send top row up
buf = u(nx+1:2*nx+1)
call MPI_Ssend(buf, nx, MPI_REAL8, up_neighbour, 0, MPI_COMM_WORLD, ierror)
! send bottom row down
buf = u(nx*ny+1:nx*ny+nx+1)
call MPI_Ssend(buf, nx, MPI_REAL8, down_neighbour, 1, MPI_COMM_WORLD, ierror)
! receive up neighbour's bottom border region
call MPI_Recv(buf, nx, MPI_REAL8, up_neighbour, 2, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierror)
u(1:nx+1) = buf
! receive down neighbour's top border region
call MPI_Recv(buf, nx, MPI_REAL8, down_neighbour, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierror)
u((ny+1)*nx+1:(ny+1)*nx+nx+1) = buf
else
! alternate order of communication for odd ranks
call MPI_Recv(buf, nx, MPI_REAL8, down_neighbour, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierror)
u((ny + 1)*nx+1:(ny + 1)*nx+nx+1) = buf
call MPI_Recv(buf, nx, MPI_REAL8, up_neighbour, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierror)
u(1:nx+1) = buf
buf = u(nx*ny+1:nx*ny+nx+1)
call MPI_Ssend(buf, nx, MPI_REAL8, down_neighbour, 2, MPI_COMM_WORLD, ierror)
buf = u(nx+1:2*nx+1)
call MPI_Ssend(buf, nx, MPI_REAL8, up_neighbour, 3, MPI_COMM_WORLD, ierror)
endif
return
end subroutine halo_exchange
subroutine evolve(u, u_prev, nx, ny)
use params
implicit none
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(inout) :: u(grid_size), u_prev(grid_size)
integer (kind=4) :: i, j, self, left, right, up, down
do i = 1,ny
do j = 0, nx-1
self = i * nx + j + 1
left = modulo(i * nx + (j - 1 + ny), ny) + 1
right = modulo(i * nx + (j + 1), ny) + 1
up = modulo((i - 1 + nx), nx) * nx + j + 1;
down = modulo((i + 1), nx) * nx + j + 1;
u(self) = u_prev(self) + dt * a * &
((u_prev(left) - 2.0 * u_prev(self) + u_prev(right)) / dx2 + &
(u_prev(up) - 2.0 * u_prev(self) + u_prev(down)) / dy2)
end do
end do
return
end subroutine evolve
subroutine iterate(u, u_prev, nx, ny, rank, unisize)
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(inout) :: u(grid_size), u_prev(grid_size)
integer (kind=4) :: i
real (kind=8) :: temp(grid_size)
do i = 1,nt
call halo_exchange(u_prev, nx, ny, rank, unisize)
call evolve(u, u_prev, nx, ny)
temp = u_prev
u_prev = u
u = temp
end do
return
end subroutine iterate
end module funcs
program heateq
use mpi
use params
use funcs
implicit none
integer :: ierror, rank, unisize
integer (kind=4) :: temp_ny, ny_local, nx_local
real (kind=8), allocatable :: init_grid(:), grid(:)
character (len=14) :: filename
call MPI_Init(ierror)
call random_seed()
call initialise()
filename = 'final_grid.txt'
call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierror)
call MPI_Comm_size(MPI_COMM_WORLD, unisize, ierror)
!
! calculate size of local domain
! distribute rows amongst process as evenly as possible
!
temp_ny = nyg / unisize
if (rank .lt. modulo(nyg,unisize)) then
temp_ny = temp_ny + 1
endif
ny_local = temp_ny
nx_local = nxg
! +2 for halo regions in y dimension
grid_size = nx_local*(ny_local+2)
allocate(init_grid(grid_size),grid(grid_size))
call randomise(init_grid)
call print_grid(init_grid, nx_local, ny_local, rank, unisize)
call iterate(grid, init_grid, nx_local, ny_local, rank, unisize)
call print_grid(grid, nx_local, ny_local, rank, unisize)
call write_grid_to_file(filename, grid, nx_local, ny_local, rank, unisize)
deallocate(init_grid,grid)
call MPI_Finalize(ierror)
end program heateq
......@@ -2,6 +2,6 @@
## Configuration
Parameters can be edited in the [param.h](./params.h) file. `nx`, `ny`, and `nt`
are the size of the grid in the x and y dimensions and the number of timesteps
respectively.
Parameters can be edited in the [param.h](./params.h) file or the main Fortran
file for C and Fortran respectively. `nx`, `ny`, and `nt` are the size of the
grid in the x and y dimensions and the number of timesteps respectively.
! initialise with randomly with even distribution
subroutine randomise(u)
use params
implicit none
real (kind=8), intent(out) :: u(grid_size)
real :: rvec(grid_size)
call random_number(rvec)
u = rvec*10.0
return
end subroutine randomise
subroutine print_local_grid(u, nx, ny)
use params
implicit none
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
integer (kind=4) :: i, j
! don't print halo regions
do i = 1,ny
do j = 0, nx-1
write(6,'(f8.3,1x)',advance='NO') u(i * nx + j+1)
end do
write(6,*)
end do
return
end subroutine print_local_grid
subroutine print_grid(u, nx, ny, rank, unisize)
use mpi
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
integer :: ierror
integer (kind=4) :: i
do i = 0,unisize - 1
if (rank .eq. i) then
call print_local_grid(u, nx, ny)
endif
call MPI_Barrier(MPI_COMM_WORLD, ierror)
end do
! last process prints an extra line
if (rank .eq. i) then
call print_local_grid(u, nx, ny)
write(6,*)
endif
call MPI_Barrier(MPI_COMM_WORLD,ierror)
return
end subroutine print_grid
subroutine write_local_grid_to_file(u, nx, ny)
use params
implicit none
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
integer (kind=4) :: i, j
do i = 1,ny
do j = 0, nx-1
write(unitnum,'(f8.3,1x)',advance='NO') u(i * nx + j+1)
end do
write(unitnum,*)
end do
return
end subroutine write_local_grid_to_file
subroutine write_grid_to_file(filename, u, nx, ny, rank, unisize)
use mpi
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(in) :: u(grid_size)
character (len=14), intent(in) :: filename
integer :: ierror
integer (kind=4) :: i
! first process truncates file to size 0 on open
if (rank .eq. 0) then
open(unit=unitnum,file=filename,access='sequential',form='FORMATTED',status='UNKNOWN',&
action='WRITE',iostat=ierror)
if (ierror .eq. 0) then
call write_local_grid_to_file(u, nx, ny)
close(unit=unitnum,status='KEEP')
endif
endif
call MPI_Barrier(MPI_COMM_WORLD, ierror)
do i = 1,unisize-1
if (rank .eq. i) then
open(unit=unitnum,file=filename,access='sequential',form='FORMATTED',status='OLD',&
action='WRITE',position='APPEND',iostat=ierror)
if (ierror .eq. 0) then
call write_local_grid_to_file(u, nx, ny)
close(unit=unitnum,status='KEEP')
endif
endif
call MPI_Barrier(MPI_COMM_WORLD, ierror)
end do
return
end subroutine write_grid_to_file
module params
integer (kind=4), parameter :: nxg = 10, nyg = 13, unitnum = 2
integer (kind=4), parameter :: nt = 1000 ! number of time-steps to evolve the system
integer (kind=4) :: grid_size
real (kind=8), parameter :: a = 1.0 ! diffusion constant
real (kind=8) :: dx, dy, dx2, dy2, dt
contains
subroutine initialise()
dx = 1.0 / real(nxg)
dy = 1.0 / real(nyg)
dx2 = dx * dx
dy2 = dy * dy
dt = dx2*dy2 / (2.0*a*(dx2 + dy2))
end subroutine initialise
end module params
module funcs
integer :: req(4)
real (kind=8), allocatable :: buf(:,:)
contains
subroutine halo_exchange(nx, ny, rank, unisize)
use mpi
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
integer :: ierror, up_neighbour, down_neighbour
up_neighbour = modulo((rank - 1 + unisize), unisize)
down_neighbour = modulo((rank + 1), unisize)
! receive up neighbour's bottom border region
call MPI_Irecv(buf(:,1), nx, MPI_REAL8, up_neighbour, 1, MPI_COMM_WORLD, req(1), ierror)
! receive down neighbour's top border region
call MPI_Irecv(buf(:,2), nx, MPI_REAL8, down_neighbour, 0, MPI_COMM_WORLD, req(2), ierror)
! send top row up
call MPI_Issend(buf(:,3), nx, MPI_REAL8, up_neighbour, 0, MPI_COMM_WORLD, req(3), ierror)
! send bottom row down
call MPI_Issend(buf(:,4), nx, MPI_REAL8, down_neighbour, 1, MPI_COMM_WORLD, req(4), ierror)
return
end subroutine halo_exchange
subroutine evolve(u, u_prev, nx, ny)
use params
implicit none
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(inout) :: u(grid_size), u_prev(grid_size)
integer (kind=4) :: i, j, self, left, right, up, down
do i = 1,ny
do j = 0, nx-1
self = i * nx + j + 1
left = modulo(i * nx + (j - 1 + ny), ny) + 1
right = modulo(i * nx + (j + 1), ny) + 1
up = modulo((i - 1 + nx), nx) * nx + j + 1;
down = modulo((i + 1), nx) * nx + j + 1;
u(self) = u_prev(self) + dt * a * &
((u_prev(left) - 2.0 * u_prev(self) + u_prev(right)) / dx2 + &
(u_prev(up) - 2.0 * u_prev(self) + u_prev(down)) / dy2)
end do
end do
return
end subroutine evolve
subroutine iterate(u, u_prev, nx, ny, rank, unisize)
use mpi
use params
implicit none
integer, intent(in) :: rank, unisize
integer (kind=4), intent(in) :: nx, ny
real (kind=8), intent(inout) :: u(grid_size), u_prev(grid_size)
integer :: ierror
integer (kind=4) :: i
real (kind=8) :: temp(grid_size)
allocate(buf(nx,4))
do i = 1,nt
buf(:,3) = u_prev(nx+1:2*nx+1)
buf(:,4) = u_prev(nx*ny+1:nx*ny+nx+1)
call halo_exchange(nx, ny, rank, unisize)
call evolve(u, u_prev, nx, ny)
call MPI_WaitAll(4,req,MPI_STATUSES_IGNORE,ierror)
u_prev(1:nxfff+1) = buf(:,1)
u((ny+1)*nx+1:(ny+1)*nx+nx+1) = buf(:,2)
temp = u_prev
u_prev = u
u = temp
end do
deallocate(buf)
return
end subroutine iterate
end module funcs
program heateq
use mpi
use params
use funcs
implicit none
integer :: ierror, rank, unisize
integer (kind=4) :: temp_ny, ny_local, nx_local
real (kind=8), allocatable :: init_grid(:), grid(:)
character (len=14) :: filename
call MPI_Init(ierror)
call random_seed()
call initialise()
filename = 'final_grid.txt'
call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierror)
call MPI_Comm_size(MPI_COMM_WORLD, unisize, ierror)
!
! calculate size of local domain
! distribute rows amongst process as evenly as possible
!
temp_ny = nyg / unisize
if (rank .lt. modulo(nyg,unisize)) then
temp_ny = temp_ny + 1
endif
ny_local = temp_ny
nx_local = nxg
! +2 for halo regions in y dimension
grid_size = nx_local*(ny_local+2)