Lisandro Dalcin Mpi4py

MPI for Python
http://mpi4py.scipy.org
Lisandro Dalcin
dalcinl@gmail.com
Centro Internacional de M
etodos Computacionales en Ingeniera
Consejo Nacional de Investigaciones Cientficas y T
ecnicas
Santa Fe, Argentina
January, 2011
Python for parallel scientific computing
PASI, Valparaso, Chile
Outline
Overview
Communicators
Point to Point Communication
Collective Operations
Compute Pi
Mandelbrot Set
Dynamic Process Management
Overview
Communicators
Compute Pi
Mandelbrot Set
What is mpi4py?
Full-featured Python bindings for MPI.

API based on the standard MPI-2 C++ bindings.
Almost all MPI calls are supported.
targeted to MPI-2 implementations.
also works with MPI-1 implementations.
Implementation
Implemented with Cython

Code base far easier to write, maintain, and extend.
Faster than other solutions (mixed Python and C codes).
A pythonic API that runs at C speed !
Features MPI-1
Process groups and communication domains.

intracommunicators
intercommunicators
Point to point communication.

blocking (send/recv)
nonblocking (isend/irecv + test/wait)
Collective operations.
Synchronization (barrier)
Communication (broadcast, scatter/gather)
Global reductions (reduce, scan)
Features MPI-2
Dynamic process management (spawn, connect/accept).

Parallel I/O (read/write).
One sided operations, a.k.a. RMA (put/get/accumulate).
Extended collective operations.
Features Python
Communication of Python objects.

high level and very convenient, based in pickle serialization
can be slow for large data (CPU and memory consuming)
<object> pickle.dump() send()
<object> pickle.load() recv()

Communication of array data (e.g. NumPy arrays).
lower level, slightly more verbose
very fast, almost C speed (for messages above 5-10 KB)
message = [<object>, (count, displ), datatype]
Point to Point Throughput Gigabit Ethernet
120
Throughput [MiB/s]
100
PingPong
Pickle
Buffer
C
80
60
40
20
0
100
101
102
103 104 105

Array Size [Bytes]
106
107
Throughput [MiB/s]
Point to Point Throughput Shared Memory
4500
4000
3500
3000
2500
2000
1500
1000
500
0
100
PingPong
Pickle
Buffer
C
101
102
103 104 105

Array Size [Bytes]
106
107
Features IPython
Integration with IPython enables MPI to be used interactively.
Start engines with MPI enabled
$ ipcluster mpiexec -n 16 --mpi=mpi4py
Connect to the engines
$ ipython
In [1]: from IPython.kernel import client
In [2]: mec = client.MultiEngineClient()
In [3]: mec.activate()
Execute commands using %px
In [4]: %px from mpi4py import MPI
In [5]: %px print(MPI.Get_processor_name())
Features Interoperability
Good support for wrapping other MPI-based codes.

You can use Cython (cimport statement).
You can use SWIG (typemaps provided).
You can use F2Py (py2f()/f2py() methods).
You can use Boost::Python or hand-written C extensions.
mpi4py will allow you to use virtually any MPI based
C/C++/Fortran code from Python.
Hello World!
from mpi4py import MPI
2
3
4
5
rank = MPI.COMM_WORLD.Get_rank()
size = MPI.COMM_WORLD.Get_size()
name = MPI.Get_processor_name()
6
7
8
9
print ("Hello, World! "

"I am process %d of %d on %s" %
(rank, size, name))
Hello World! Wrapping with SWIG

C source
1
2
3
4
5
6
7
8
9
10
11
/* file: helloworld.c */
void sayhello(MPI_Comm comm)
{
int size, rank;
MPI_Comm_size(comm, &size);
MPI_Comm_rank(comm, &rank);
printf("Hello, World! "
"I am process "
"%d of %d.\n",
rank, size);
}
SWIG interface file

1
2
3
4
5
6
// file: helloworld.i
%module helloworld
%{
# include <mpi.h>
# include "helloworld.c"
}%
7
8
9
%include mpi4py/mpi4py.i
%mpi4py_typemap(Comm, MPI_Comm);
10
11
void sayhello(MPI_Comm comm);
At the Python prompt . . .

>>> from mpi4py import MPI
>>> import helloworld
>>> helloworld.sayhello(MPI.COMM_WORLD)
Hello, World! I am process 0 of 1.
Hello World! Wrapping with Boost.Python

1
2
3
4
// file: helloworld.cxx
# include <boost / python.hpp>
# include <mpi4py / mpi4py.h>
using namespace boost::python;
5
6
7
8
9
10
11
12
# include "helloworld.c"
static void wrap_sayhello(object py_comm) {
PyObject* py_obj = py_comm.ptr();
MPI_Comm *comm_p = PyMPIComm_Get(py_obj);
if (comm_p == NULL) throw_error_already_set();
sayhello(*comm_p);
}
13
14
15
16
17
BOOST_PYTHON_MODULE(helloworld) {
if (import_mpi4py() < 0) return;
def("sayhello", wrap_sayhello);
}
Hello World! Wrapping with F2Py

Fortran 90 source
1
2
3
4
5
6
7
8
9
! file: helloworld.f90
subroutine sayhello(comm)
use mpi
implicit none
integer :: comm, rank, size, ierr
call MPI_Comm_size(comm, size, ierr)
call MPI_Comm_rank(comm, rank, ierr)
print *, Hello, World! I am process ,rank, of ,size,.
end subroutine sayhello
At the Python prompt . . .

>>> from mpi4py import MPI
>>> import helloworld
>>> fcomm = MPI.COMM_WORLD.py2f()
>>> helloworld.sayhello(fcomm)
Hello, World! I am process 0 of 1.
Overview
Communicators
Compute Pi
Mandelbrot Set
Communicators
communicator = process group + communication context
Predefined instances
COMM WORLD
COMM SELF
COMM NULL
Accessors
rank = comm.Get rank() # or comm.rank
size = comm.Get size() # or comm.size
group = comm.Get group()
Constructors
newcomm = comm.Dup()
newcomm = comm.Create(group)
newcomm = comm.Split(color, key)
Communicators Create()
1
2
3
4
comm = MPI.COMM_WORLD
group = comm.Get_group()
5
6
7
newgroup = group.Excl([0])
newcomm = comm.Create(newgroup)
8
9
10
11
12
13
if comm.rank == 0:
assert newcomm == MPI.COMM_NULL
else:
assert newcomm.size == comm.size - 1
assert newcomm.rank == comm.rank - 1
14
15
16
group.Free(); newgroup.Free()
if newcomm: newcomm.Free()
Communicators Split()
1
2
3
4
world_rank = MPI.COMM_WORLD.Get_rank()
world_size = MPI.COMM_WORLD.Get_size()
5
6
7
8
9
10
11
if world_rank < world_size//2:

color = 55
key = -world_rank
else:
color = 77
key = +world_rank
12
13
14
15
newcomm = MPI.COMM_WORLD.Split(color, key)

# ...
newcomm.Free()
Exercise #1
a) Create a new process group containing the processes in the

group of COMM WORLD with even rank. Use the new group to
create a new communicator.
Tip: use Group.Incl() or Group.Range incl()
b) Use Comm.Split() to split COMM WORLD in two halves.
The first half contains the processes with even rank in
COMM WORLD. The process rank ordering in the new
communication is ascending.
The second half contains the processes with odd rank in
COMM WORLD. The process rank ordering in the new
communication is descending.
Overview
Communicators
Compute Pi
Mandelbrot Set
Blocking communication
Python objects
comm.send(obj, dest=0, tag=0)
obj = comm.recv(None, src=0, tag=0)
Array data
comm.Send([array, count, datatype], dest=0, tag=0)
comm.Recv([array, count, datatype], src=0, tag=0)
Nonblocking communication
Python objects
request = comm.isend(object, dest=0, tag=0)}
request.Wait()
Array data
req1 = comm.Isend([array, count, datatype], dest=0, tag=0)
req2 = comm.Irecv([array, count, datatype], src=0, tag=0)
MPI.Request.Waitall([req1, req2])}
PingPong
1
2
3

assert comm.size == 2
4
5
6
7
8
9
10
11
12
if comm.rank == 0:
sendmsg = 777
comm.send(sendmsg, dest=1, tag=55)
recvmsg = comm.recv(source=1, tag=77)
else:
recvmsg = comm.recv(source=0, tag=55)
sendmsg = "abc"
comm.send(sendmsg, dest=0, tag=77)
PingPing
1
2
3

4
5
6
7
8
9
10
if comm.rank == 0:
sendmsg = 777
target = 1
else:
sendmsg = "abc"
target = 0
11
12
13
14
request = comm.isend(sendmsg, dest=target, tag=77)

recvmsg = comm.recv(source=target, tag=77)
request.Wait()
Exchange
1
2

3
4
5
6
sendmsg = [comm.rank]*3
right = (comm.rank + 1) % comm.size
left = (comm.rank - 1) % comm.size
7
8
9
10
11
req1
req2
lmsg
rmsg
=
=
=
=
comm.isend(sendmsg, dest=right)
comm.isend(sendmsg, dest=left)
comm.recv(source=left)
comm.recv(source=right)
12
13
14
15
MPI.Request.Waitall([req1, req2])
assert lmsg == [left] * 3
assert rmsg == [right] * 3
PingPing with NumPy arrays

1
2
3
4

import numpy
5
6
7
8
9
10
11
12
13
if comm.rank
array1 =
array2 =
target =
else:
array1 =
array2 =
target =
== 0:
numpy.arange(10000, dtype=f)
numpy.empty(10000, dtype=f)
1
numpy.ones(10000, dtype=f)
numpy.empty(10000, dtype=f)
0
14
15
16
17
request = comm.Isend([array1, MPI.FLOAT], dest=target)

comm.Recv([array2, MPI.FLOAT], source=target)
request.Wait()
Exercise #2
a) Modify PingPong example to communicate NumPy arrays.

Tip: use Comm.Send() and Comm.Recv()
b) Modify Exchange example to communicate NumPy arrays.
Use nonblocking communication for both sending and receiving.
Tip: use Comm.Isend() and Comm.Irecv()
Overview
Communicators
Compute Pi
Mandelbrot Set
processes
data
A
broadcast
A
A
A
A
A
0
0
0
0
0
0
data
processes
scatter
A
A
A
gather
A
A
A
0
1
2
3
4
5
processes
data
A
B
C
D
E
F
allgather
A
A
0
0
0
0
0
0
B
B
B
B
B
B
0
0
0
0
0
0
C
C
C
C
C
C
0
0
0
0
0
0
D
D
D
D
D
D
0
0
0
0
0
0
E
E
E
E
E
E
0
0
0
0
0
0
F
F
F
F
F
F
0
0
0
0
0
0
data
processes
A
B
C
D
E
F
0
0
0
0
0
0
A
B
C
D
E
F
1
1
1
1
1
1
A
B
C
D
E
F
2
2
2
2
2
2
A
B
C
D
E
F
3
3
3
3
3
3
A
B
C
D
E
F
4
4
4
4
4
4
A
B
C
D
E
F
5
5
5
5
5
5
alltoall
A
A
A
A
A
0
1
2
3
4
5
B
B
B
B
B
B
0
1
2
3
4
5
C
C
C
C
C
C
0
1
2
3
4
5
D
D
D
D
D
D
0
1
2
3
4
5
E
E
E
E
E
E
0
1
2
3
4
5
F
F
F
F
F
F
0
1
2
3
4
5
Broadcast
1
2

3
4
5
6
7
if comm.rank == 0:
sendmsg = (7, "abc", [1.0,2+3j], {3:4})
else:
sendmsg = None
8
9
recvmsg = comm.bcast(sendmsg, root=0)
Scatter
1
2

3
4
5
6
7
if comm.rank == 0:
sendmsg = [i**2 for i in range(comm.size)]
else:
sendmsg = None
8
9
recvmsg = comm.scatter(sendmsg, root=0)
Gather & Gather to All
1
2

3
4
sendmsg = comm.rank**2
5
6
recvmsg1 = comm.gather(sendmsg, root=0)
7
8
recvmsg2 = comm.allgather(sendmsg)
Reduce & Reduce to All
1
2

3
4
sendmsg = comm.rank
5
6
recvmsg1 = comm.reduce(sendmsg, op=MPI.SUM, root=0)
7
8
recvmsg2 = comm.allreduce(sendmsg)
Exercise #3
a) Modify Broadcast, Scatter, and Gather example to

communicate NumPy arrays.
b) Write a routine implementing parallel matrixvector product
y = matvec(comm,A,x).
the global matrix is dense and square.
matrix rows and vector entries have matching block distribution.
all processes own the same number of matrix rows.
Tip: use Comm.Allgather() and numpy.dot()
Overview
Communicators
Compute Pi
Mandelbrot Set
Compute Pi
=
0
1
4
dx
2
1+x
n
n1
i=0
4
2
1 + ( i+0.5
n )
Compute Pi sequential
1
import math
2
3
4
5
6
7
8
9
def compute_pi(n):
h = 1.0 / n
s = 0.0
for i in range(n):
x = h * (i + 0.5)
s += 4.0 / (1.0 + x**2)
return s * h
10
11
12
13
n = 10
pi = compute_pi(n)
error = abs(pi - math.pi)
14
15
16
print ("pi is approximately %.16f, "

"error is %.16f" % (pi, error))
Compute Pi parallel [1]

1
2

import math
3
4
5
6
7
8
9
10
def compute_pi(n, start=0, step=1):

h = 1.0 / n
s = 0.0
for i in range(start, n, step):
x = h * (i + 0.5)
s += 4.0 / (1.0 + x**2)
return s * h
11
12
13
14
nprocs = comm.Get_size()
myrank = comm.Get_rank()
Compute Pi parallel [2]

1
2
3
4
if myrank == 0:
n = 10
else:
n = None
5
6
n = comm.bcast(n, root=0)
7
8
mypi = compute_pi(n, myrank, nprocs)
9
10
pi = comm.reduce(mypi, op=MPI.SUM, root=0)
11
12
13
14
15
if myrank == 0:
error = abs(pi - math.pi)
"error is %.16f" % (pi, error))
Exercise #4
Modify Compute Pi example to employ NumPy.

Tip: you can convert a Python int/float object to a NumPy
scalar with x = numpy.array(x).
Overview
Communicators
Compute Pi
Mandelbrot Set
Mandelbrot Set
Mandelbrot Set sequential [1]

1
2
3
4
5
6
7
8
def mandelbrot(x, y, maxit):

c = x + y*1j
z = 0 + 0j
it = 0
while abs(z) < 2 and it < maxit:
z = z**2 + c
it += 1
return it
9
10
11
12
13
x1, x2
y1, y2
w, h
maxit
=
=
=
=
-2.0, 1.0
-1.0, 1.0
150, 100
127
Mandelbrot Set sequential [2]

1
2
3
4
5
6
7
8
9
import numpy
C = numpy.zeros([h, w], dtype=i)
dx = (x2 - x1) / w
dy = (y2 - y1) / h
for i in range(h):
y = y1 + i * dy
for j in range(w):
x = x1 + j * dx
C[i, j] = mandelbrot(x, y, maxit)
10
11
12
13
14
from matplotlib import pyplot

pyplot.imshow(C, aspect=equal)
pyplot.spectral()
pyplot.show()
Mandelbrot Set partitioning
0
1
2
Block distribution
0
0
0
2
Cyclic distribution
Mandelbrot Set parallel, block [1]

1
2
3
4
5
6
7
8
def mandelbrot(x, y, maxit):

c = x + y*1j
z = 0 + 0j
it = 0
while abs(z) < 2 and it < maxit:
z = z**2 + c
it += 1
return it
9
10
11
12
13
x1, x2
y1, y2
w, h
maxit
=
=
=
=
-2.0, 1.0
-1.0, 1.0
150, 100
127

1
2

import numpy
3
4
5
6
size = comm.Get_size()
rank = comm.Get_rank()
7
8
9
# number of rows to compute here

N = h // size + (h % size > rank)
10
11
12
# first row to compute here

start = comm.scan(N)-N
13
14
15
# array to store local result

Cl = numpy.zeros([N, w], dtype=i)

1
# compute owned rows
2
3
4
5
6
7
8
9
dx = (x2 - x1) / w
dy = (y2 - y1) / h
for i in range(N):
y = y1 + (i + start) * dy
for j in range(w):
x = x1 + j * dx
Cl[i, j] = mandelbrot(x, y, maxit)

1
# gather results at root (process 0)
2
3
4
5
6
counts = comm.gather(N, root=0)

C = None
if rank == 0:
C = numpy.zeros([h, w], dtype=i)
7
8
9
rowtype = MPI.INT.Create_contiguous(w)
rowtype.Commit()
10
11
12
13
comm.Gatherv(sendbuf=[Cl, MPI.INT],
recvbuf=[C, (counts, None), rowtype],
root=0)
14
15
rowtype.Free()

1
2
3
4
5
if comm.rank == 0:
from matplotlib import pyplot
pyplot.imshow(C, aspect=equal)
pyplot.spectral()
pyplot.show()
Exercise #5
Measure the wall clock time Ti of local computations at each

process for the Mandelbrot Set example with block and cyclic row
distributions.
What row distribution is better regarding load balancing?
Tip: use Wtime() to measure wall time, compute the ratio
Tmax /Tmin to compare load balancing.
Overview
Communicators
Compute Pi
Mandelbrot Set
Useful in assembling complex distributed applications. Can

couple independent parallel codes written in different
languages.
Create new processes from a running program.
Comm.Spawn() and Comm.Get parent()
Connect two running applications together.
Comm.Connect() and Comm.Accept()
Dynamic Process Management Spawning

Spawning new processes is a collective operation that creates an
intercommunicator.
Local group is group of spawning processes (parent).
Remote group is group of new processes (child).
Comm.Spawn() lets parent processes spawn the child
processes. It returns a new intercommunicator.
Comm.Get parent() lets child processes find
intercommunicator to the parent group. Child processes have
own COMM WORLD.
Comm.Disconnect() ends the parentchild connection. After
that, both groups can continue running.
Dynamic Process Management Compute Pi (parent)

1
2

import sys, numpy
3
4
5
6
comm = MPI.COMM_SELF.Spawn(sys.executable,
args=[compute_pi-child.py],
maxprocs=5)
7
8
9
10
11
12
13
N = numpy.array(10, i)
comm.Bcast([N, MPI.INT], root=MPI.ROOT)
PI = numpy.array(0.0, d)
comm.Reduce(None, [PI, MPI.DOUBLE],
op=MPI.SUM, root=MPI.ROOT)
comm.Disconnect()
14
15
16
17
error = abs(PI - numpy.pi)

"error is %.16f" % (PI, error))
Dynamic Process Management Compute Pi (child)

1
2

import numpy
3
4
5
6
comm = MPI.Comm.Get_parent()
size = comm.Get_size()
rank = comm.Get_rank()
7
8
9
10
11
12
13
14
15
16
N = numpy.array(0, dtype=i)
comm.Bcast([N, MPI.INT], root=0)
h = 1.0 / N; s = 0.0
for i in range(rank, N, size):
x = h * (i + 0.5)
s += 4.0 / (1.0 + x**2)
PI = numpy.array(s * h, dtype=d)
comm.Reduce([PI, MPI.DOUBLE], None,
op=MPI.SUM, root=0)
17
18
comm.Disconnect()
Exercise #5
a) Implement the Compute Pi child code in C or C++ . Adjust

the parent code in Python to spawn the new implementation.
b) Compute and plot the Mandelbrot Set using spawning with
parent/child codes implemented in Python.
Tip: Reuse the provided parent code in Python and translate
the child code in Fortran 90 to Python.
Do not hesitate to ask for help . . .

Mailing List: mpi4py@googlegroups.com
Mail&Chat: dalcinl@gmail.com
Thanks!

Lisandro Dalcin Mpi4py

Uploaded by

Copyright:

Available Formats

Lisandro Dalcin Mpi4py

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Lisandro Dalcin Mpi4py

Uploaded by

Copyright:

Available Formats

MPI for Python

Full-featured Python bindings for MPI.

Implemented with Cython

Process groups and communication domains.

Point to point communication.

Dynamic process management (spawn, connect/accept).

Communication of Python objects.

<object> pickle.dump() send()

<object> pickle.load() recv()

message = [<object>, (count, displ), datatype]

Point to Point Throughput Gigabit Ethernet

103 104 105

Point to Point Throughput Shared Memory

103 104 105

Good support for wrapping other MPI-based codes.

from mpi4py import MPI

print ("Hello, World! "

Hello World! Wrapping with SWIG

SWIG interface file

void sayhello(MPI_Comm comm);

At the Python prompt . . .

Hello World! Wrapping with Boost.Python

Hello World! Wrapping with F2Py

At the Python prompt . . .

from mpi4py import MPI

from mpi4py import MPI

if world_rank < world_size//2:

newcomm = MPI.COMM_WORLD.Split(color, key)

a) Create a new process group containing the processes in the

from mpi4py import MPI

from mpi4py import MPI

request = comm.isend(sendmsg, dest=target, tag=77)

from mpi4py import MPI

PingPing with NumPy arrays

from mpi4py import MPI

request = comm.Isend([array1, MPI.FLOAT], dest=target)

a) Modify PingPong example to communicate NumPy arrays.

from mpi4py import MPI

recvmsg = comm.bcast(sendmsg, root=0)

from mpi4py import MPI

recvmsg = comm.scatter(sendmsg, root=0)

Gather & Gather to All

from mpi4py import MPI

recvmsg1 = comm.gather(sendmsg, root=0)

Reduce & Reduce to All

from mpi4py import MPI

recvmsg1 = comm.reduce(sendmsg, op=MPI.SUM, root=0)

a) Modify Broadcast, Scatter, and Gather example to

Tip: use Comm.Allgather() and numpy.dot()

print ("pi is approximately %.16f, "

Compute Pi parallel [1]

from mpi4py import MPI

def compute_pi(n, start=0, step=1):

Compute Pi parallel [2]

mypi = compute_pi(n, myrank, nprocs)

pi = comm.reduce(mypi, op=MPI.SUM, root=0)

Modify Compute Pi example to employ NumPy.

Mandelbrot Set sequential [1]

def mandelbrot(x, y, maxit):