Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Programming with MPI

OpenMP lets you use all the cores of a single node, but the largest scientific simulations run on thousands of nodes. To go beyond one node, processes have to communicate explicitly by sending messages to each other over the network. The Message Passing Interface (MPI) is the standard library for doing this.

In this lecture, we will:

Parallel programming languages

Evolution of programming methods:

Distributed memory versus shared memory

MPI history

MPI implementations

Open source libraries can be installed on almost any architecture (for example on your laptop):

Vendors also provide their own implementations, for example:

MPI tools

Debuggers and performance analysis tools:

Scientific libraries:

General concepts

Parallel processing

MPI is a library which allows process coordination and scheduling between millions of processors using a message-passing paradigm.

Message attributes

The MPI environment

Data distribution

Data (grid cells or particles) are distributed between nodes and cores using a domain decomposition strategy. On the left, a regular grid is split into 9 equal domains. On the right, the domains of a particle simulation adapt to the particle density, so that each process gets roughly the same amount of work.

MPI basics

Setting up the MPI environment

The MPI definitions are made available with:

The MPI environment is launched with the MPI_INIT() routine and terminated with the MPI_FINALIZE() routine.

C syntax:

int MPI_Init(int *argc, char ***argv);
int MPI_Finalize(void);

Fortran syntax:

call MPI_INIT(code)
call MPI_FINALIZE(code)

Who am I?

Each process can query the total number of processes in the communicator MPI_COMM_WORLD with MPI_COMM_SIZE(), and its own identifier (its rank) with MPI_COMM_RANK():

program who_am_I
  use mpi
  implicit none
  integer :: nb_procs,rank,code

  call MPI_INIT(code)

  call MPI_COMM_SIZE(MPI_COMM_WORLD,nb_procs,code)
  call MPI_COMM_RANK(MPI_COMM_WORLD,rank,code)

  print *,'I am the process ',rank,' among ',nb_procs

  call MPI_FINALIZE(code)
end program who_am_I
$ mpif90 who_am_I.f90 -o who_am_I
$ mpiexec -n 7 who_am_I
 I am the process            3  among            7
 I am the process            0  among            7
 I am the process            4  among            7
 I am the process            1  among            7
 I am the process            5  among            7
 I am the process            2  among            7
 I am the process            6  among            7

The processes print in a random order: they all run independently.

Point-to-point communications

Blocking send and receive

program point_to_point
  use mpi
  implicit none

  integer, dimension(MPI_STATUS_SIZE) :: status
  integer, parameter                  :: tag=100
  integer                             :: rank,value,code

  call MPI_INIT(code)

  call MPI_COMM_RANK(MPI_COMM_WORLD,rank,code)

  if (rank == 2) then
     value=1000
     call MPI_SEND(value,1,MPI_INTEGER,5,tag,MPI_COMM_WORLD,code)
  elseif (rank == 5) then
     call MPI_RECV(value,1,MPI_INTEGER,2,tag,MPI_COMM_WORLD,status,code)
     print *,'I, process 5, I received ',value,' from the process 2'
  end if

  call MPI_FINALIZE(code)

end program point_to_point
$ mpiexec -n 7 point_to_point
 I, process 5, I received         1000  from the process 2

With blocking send and receive, if you are not careful, you can have a deadlock. This means the processes will wait for something that will never happen. The classical mistake is to first call a blocking send towards the next process, and then call a blocking receive from the previous process.

Non-blocking send and receive

Communication costs can be large. The InfiniBand network latency is around a microsecond, equivalent to thousands of processor cycles. One has to add the cost of the message itself (limited by the bandwidth of the InfiniBand network, around 10 GB/s).

Non-blocking routines can be used to perform calculations in parallel to the communication (an additional level of parallelism). There is no risk of deadlock, but there is a risk of memory leak if the communication is not properly terminated.

The communication buffers, however, cannot be used before the proper completion of the send or receive. The user needs to test that the communication is successful in order to proceed with the data, which leads to a higher algorithmic complexity.

Computation-communication overlap.

Example: halo exchange

A classical use case is a stencil computation on a 2D grid decomposed into domains. Each process needs the values of the cells along the boundaries of its 4 neighbors (North, South, West, and East), which are stored in extra “ghost” cells around its own domain.

The communications with the 4 neighbors are started with non-blocking calls. The derived datatypes rowtype and columntype describe a row and a column of the array, and the array indices of the first element of each row or column are omitted here for brevity:

SUBROUTINE start_communication(u)
  ! Send to the North and receive from the South
  CALL MPI_IRECV( u(,), 1, rowtype, neighbor(S), &
       tag, comm2d, request(1), code)
  CALL MPI_ISEND( u(,), 1, rowtype, neighbor(N), &
       tag, comm2d, request(2), code)

  ! Send to the South and receive from the North
  CALL MPI_IRECV( u(,), 1, rowtype, neighbor(N), &
       tag, comm2d, request(3), code)
  CALL MPI_ISEND( u(,), 1, rowtype, neighbor(S), &
       tag, comm2d, request(4), code)

  ! Send to the West and receive from the East
  CALL MPI_IRECV( u(,), 1, columntype, neighbor(E), &
       tag, comm2d, request(5), code)
  CALL MPI_ISEND( u(,), 1, columntype, neighbor(W), &
       tag, comm2d, request(6), code)

  ! Send to the East and receive from the West
  CALL MPI_IRECV( u(,), 1, columntype, neighbor(W), &
       tag, comm2d, request(7), code)
  CALL MPI_ISEND( u(,), 1, columntype, neighbor(E), &
       tag, comm2d, request(8), code)
END SUBROUTINE start_communication

SUBROUTINE end_communication(u)
  CALL MPI_WAITALL(2*NB_NEIGHBORS, request, tab_status, code)
END SUBROUTINE end_communication

While the messages are in flight, each process updates the interior of its domain, which does not depend on the ghost cells. Once the communications are completed, it updates the 4 boundaries:

DO WHILE ((.NOT. convergence) .AND. (it < it_max))
  it = it +1
  u(sx:ex,sy:ey) = u_new(sx:ex,sy:ey)

  ! Exchange value on the interfaces
  CALL start_communication( u )

  ! Compute u
  CALL calcul( u, u_new, sx+1, ex-1, sy+1, ey-1)

  CALL end_communication( u )

  ! North
  CALL calcul( u, u_new, sx, sx, sy, ey)
  ! South
  CALL calcul( u, u_new, ex, ex, sy, ey)
  ! West
  CALL calcul( u, u_new, sx, ex, sy, sy)
  ! East
  CALL calcul( u, u_new, sx, ex, ey, ey)

  ! Compute global error
  diffnorm = global_error (u, u_new)

  convergence = ( diffnorm < eps )

END DO

Collective communications

General concepts

Types of collectives

Scatter

Process 2 splits its array of 8 values into blocks of 2 values, and sends one block to each process (including itself):

program scatter
  use mpi
  implicit none

  integer, parameter                :: nb_values=8
  integer                           :: nb_procs,rank,block_length,i,code
  real, allocatable, dimension(:)   :: values,data

  call MPI_INIT(code)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,nb_procs,code)
  call MPI_COMM_RANK(MPI_COMM_WORLD,rank,code)
  block_length=nb_values/nb_procs
  allocate(data(block_length))

  if (rank == 2) then
     allocate(values(nb_values))
     values(:)=(/(1000.+i,i=1,nb_values)/)
     print *,'I, process ',rank,'send my values array : ',&
             values(1:nb_values)
  end if

  call MPI_SCATTER(values,block_length,MPI_REAL,data,block_length, &
                   MPI_REAL,2,MPI_COMM_WORLD,code)
  print *,'I, process ',rank,', received ', data(1:block_length), &
          ' of process 2'
  call MPI_FINALIZE(code)

end program scatter
$ mpiexec -n 4 scatter
 I, process            2 send my values array :    1001.00000       1002.00000       1003.00000       1004.00000       1005.00000       1006.00000       1007.00000       1008.00000
 I, process            0 , received    1001.00000       1002.00000      of process 2
 I, process            1 , received    1003.00000       1004.00000      of process 2
 I, process            3 , received    1007.00000       1008.00000      of process 2
 I, process            2 , received    1005.00000       1006.00000      of process 2

All reduce

Each process contributes a value, and all processes receive the product of all the values, here 10×1×2×3×4×5×6=720010 \times 1 \times 2 \times 3 \times 4 \times 5 \times 6 = 7200:

program allreduce

  use mpi
  implicit none

  integer :: nb_procs,rank,value,product,code

  call MPI_INIT(code)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,nb_procs,code)
  call MPI_COMM_RANK(MPI_COMM_WORLD,rank,code)

  if (rank == 0) then
     value=10
  else
     value=rank
  endif

  call MPI_ALLREDUCE(value,product,1,MPI_INTEGER,MPI_PROD,MPI_COMM_WORLD,code)

  print *,'I, process ',rank,', received the value of the global product ',product

  call MPI_FINALIZE(code)

end program allreduce
$ mpiexec -n 7 allreduce
 I, process            6 , received the value of the global product         7200
 I, process            2 , received the value of the global product         7200
 I, process            0 , received the value of the global product         7200
 I, process            4 , received the value of the global product         7200
 I, process            5 , received the value of the global product         7200
 I, process            3 , received the value of the global product         7200
 I, process            1 , received the value of the global product         7200

OpenMP versus MPI

Graphics Processing Units (GPU)

GPU model

A CPU core runs 1 thread at a time, with fast access to its main memory. A GPU runs about 32 threads per core, with very fast access to its own memory. But data has to be copied from the CPU main memory to the GPU memory and back, which is slow.

GPU programming

Conclusion

Parallel computing methods:

Week 7: Parallel Computing
Parallel concepts
WIP: Week 8: GPU Programming
Programming GPUs with CUDA