{"id":471,"date":"2019-08-02T14:27:52","date_gmt":"2019-08-02T12:27:52","guid":{"rendered":"https:\/\/www.kth.se\/blogs\/pdc\/?p=471"},"modified":"2019-08-02T14:27:52","modified_gmt":"2019-08-02T12:27:52","slug":"parallel-programming-in-python-mpi4py-part-1","status":"publish","type":"post","link":"https:\/\/www.kth.se\/blogs\/pdc\/2019\/08\/parallel-programming-in-python-mpi4py-part-1\/","title":{"rendered":"Parallel programming in Python: mpi4py (part 1)"},"content":{"rendered":"

In previous posts<\/a> we have introduced the multiprocessing<\/code> module which makes it possible to parallelize Python programs on shared memory systems. The limitation of the multiprocessing<\/code> module is that it does not support parallelization over multiple compute nodes (i.e. on distributed memory systems). To overcome this limitation and enable cross-node parallelization, we can use MPI for Python<\/a>, that is, the mpi4py<\/code> module. This module provides an object-oriented interface that resembles the message passing interface (MPI)<\/a>,\u00a0 and hence allows Python programs to exploit multiple processors on multiple compute nodes. The mpi4py<\/code><\/a> module supports both point-to-point and collective communications for Python objects as well as buffer-like objects. This post will briefly introduce the use of the mpi4py<\/code> module in communicating generic Python objects, via all-lowercase methods including send<\/code>, recv<\/code>, isend<\/code>, irecv<\/code>, bcast<\/code>, scatter<\/code>, gather<\/code>, and reduce<\/code>.<\/p>\n

<\/p>\n

Basics of mpi4py<\/h2>\n

The mpi4py<\/code> module has been installed on Beskow<\/a>. To use mpi4py<\/code> on Beskow, we need to load the module first<\/p>\n

\n
\n
# on Beskow\r\nmodule load mpi4py\/3.0.2\/py37<\/pre>\n<\/div>\n<\/div>\n
The mpi4py\/3.0.2\/py37<\/code> module will automatically load the anaconda\/2019.03\/py37<\/code> module, which is an open-source distribution of Python for scientific computing (see Anaconda<\/a>). We recommend Python 3 because Python 2 is going to be retired soon<\/a>. After loading the mpi4py<\/code> module, we can check that the python3<\/code> command (which is needed to run\u00a0 mpi4py<\/code> code) is available in the PATH<\/code> environment variable, using the which<\/code> command<\/p>\n
\n
\n
$ which python3\r\n\/pdc\/vol\/anaconda\/2019.03\/py37\/bin\/python3<\/pre>\n<\/div>\n<\/div>\n
Now we can write our first Python program with mpi4py<\/code><\/p>\n
\n
\n
from<\/span> mpi4py import<\/span> MPI\r\n\r\ncomm =<\/span> MPI.<\/span>COMM_WORLD\r\nrank =<\/span> comm.<\/span>Get_rank(<\/span>)<\/span>\r\nsize =<\/span> comm.<\/span>Get_size(<\/span>)<\/span>\r\n\r\nprint<\/span>(<\/span>'Hello from process {} out of {}'<\/span>.<\/span>format(<\/span>rank,<\/span> size)<\/span>)<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n
To understand the code we need to know several concepts in MPI. In parallel programming with MPI, we need the so-called communicator<\/strong>, which is a group of processes that can talk to each other. To identify the processes with that group, each process is assigned a rank<\/strong> that is unique within the communicator. It also makes sense to know the total number of processes, which is often referred to as the size<\/strong> of the communicator.<\/p>\n
In the above code, we defined the variable comm<\/code> as the default communicator MPI.COMM_WORLD<\/code>. You may recognize that MPI.COMM_WORLD<\/code> in mpi4py<\/code> corresponds to MPI_COMM_WORLD<\/code> in MPI programs written in C<\/a>. The rank of each process is retrieved via the Get_rank<\/code> method of the communicator, and the size of the communicator is obtained by the Get_size<\/code> method. Now that the rank and size are known, we can print a \u201cHello\u201d message from each process.<\/p>\n
The usual way to execute an mpi4py<\/code> code in parallel is to use mpirun<\/code> and python3<\/code>, for example \u201cmpirun -n 4 python3 hello.py<\/code>\u201d will run the code on 4 processes, assuming that the code is saved in a file named \u201chello.py\u201d. On Beskow<\/a>, however, the setup is different since the resources (compute nodes) are managed by the SLURM workload manager<\/a>. We’ll need to first request one or more compute nodes using salloc<\/code>, and then execute \u201cpython3 hello.py<\/code>\u201d via srun<\/code>. For example, the output from running the code on 4 processes is as follows:<\/p>\n
\n
\n
$ salloc --nodes 1 -t 00:10:00 -A <your-project-account>\r\nsalloc: Granted job allocation <your-job-id>\r\n\r\n$ srun -n 4 python3 hello.py\r\nHello from process 0 out of 4\r\nHello from process 1 out of 4\r\nHello from process 2 out of 4\r\nHello from process 3 out of 4<\/pre>\n<\/div>\n<\/div>\n
Point-to-point communication<\/h2>\n
For point-to-point communication between Python objects, mpi4py<\/code> provides the send<\/code>\u00a0and recv<\/code> methods that are similar to those in MPI. An example of code for passing (which is usually referred to as \u201ccommunicating\u201d) a Python dictionary object between the master process (which has rank = 0) and the worker processes (that all have rank > 0) is given below.<\/p>\n
\n
\n
from<\/span> mpi4py import<\/span> MPI\r\n\r\ncomm =<\/span> MPI.<\/span>COMM_WORLD\r\nrank =<\/span> comm.<\/span>Get_rank(<\/span>)<\/span>\r\nsize =<\/span> comm.<\/span>Get_size(<\/span>)<\/span>\r\n\r\n# master process<\/span>\r\nif<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    data =<\/span> {<\/span>'x'<\/span>:<\/span> 1<\/span>,<\/span> 'y'<\/span>:<\/span> 2.0<\/span>}<\/span>\r\n    # master process sends data to worker processes by<\/span>\r\n    # going through the ranks of all worker processes<\/span>\r\n    for<\/span> i in<\/span> range<\/span>(<\/span>1<\/span>,<\/span> size)<\/span>:<\/span>\r\n        comm.<\/span>send(<\/span>data,<\/span> dest=<\/span>i,<\/span> tag=<\/span>i)<\/span>\r\n        print<\/span>(<\/span>'Process {} sent data:'<\/span>.<\/span>format(<\/span>rank)<\/span>,<\/span> data)<\/span>\r\n\r\n# worker processes<\/span>\r\nelse<\/span>:<\/span>\r\n    # each worker process receives data from master process<\/span>\r\n    data =<\/span> comm.<\/span>recv(<\/span>source=<\/span>0<\/span>,<\/span> tag=<\/span>rank)<\/span>\r\n    print<\/span>(<\/span>'Process {} received data:'<\/span>.<\/span>format(<\/span>rank)<\/span>,<\/span> data)<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n
In this example (which we will name \u201csend-recv.py\u201d), the dictionary object data<\/code> is sent from the master process to each worker process. We check the value of rank<\/code> to determine if a process is the master process. The for<\/code> loop is only executed for the master process, and i<\/code> (starting from 1) goes through the ranks of the worker processes. We can also use tag<\/code> in the send<\/code>\u00a0and recv<\/code> methods to distinguish between the messages if there are multiple communications between two processes. The output from running the example on 4 processes is as follows:<\/p>\n
\n
\n
$ srun -n 4 python3 send-recv.py\r\nProcess 0 sent data: {'x': 1, 'y': 2.0}\r\nProcess 0 sent data: {'x': 1, 'y': 2.0}\r\nProcess 0 sent data: {'x': 1, 'y': 2.0}\r\nProcess 1 received data: {'x': 1, 'y': 2.0}\r\nProcess 2 received data: {'x': 1, 'y': 2.0}\r\nProcess 3 received data: {'x': 1, 'y': 2.0}<\/pre>\n<\/div>\n<\/div>\n
Note, however, that the above example uses blocking communication methods, which means that the execution of code will not proceed until the communication is completed. This blocking behaviour is not always desirable in parallel programming; in some cases it is beneficial to use non-blocking communication methods. As shown in the code below, isend<\/code> and irecv<\/code> are non-blocking methods that immediately return Request<\/code> objects, and we can use the\u00a0wait<\/code> method to manage the completion of the communication. If necessary, one can do some computations between comm.isend(...)<\/code> and req.wait()<\/code> to increase the efficiency by overlapping the communications and computations.<\/p>\n
\n
\n
if<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    data =\u00a0{<\/span>'x'<\/span>: 1<\/span>, 'y'<\/span>: 2.0<\/span>}<\/span><\/span>\r\n    for<\/span> i in<\/span> range<\/span>(<\/span>1<\/span>,<\/span> size)<\/span>:<\/span>\r\n        req =<\/span> comm.<\/span>isend(<\/span>data,<\/span> dest=<\/span>i,<\/span> tag=<\/span>i)<\/span>\r\n        req.<\/span>wait(<\/span>)<\/span>\r\n        print<\/span>(<\/span>'Process {} sent data:'<\/span>.<\/span>format(<\/span>rank)<\/span>,<\/span> data)<\/span>\r\n\r\nelse<\/span>:<\/span>\r\n    req =<\/span> comm.<\/span>irecv(<\/span>source=<\/span>0<\/span>,<\/span> tag=<\/span>rank)<\/span>\r\n    data =<\/span> req.<\/span>wait(<\/span>)<\/span>\r\n    print<\/span>(<\/span>'Process {} received data:'<\/span>.<\/span>format(<\/span>rank)<\/span>,<\/span> data)<\/span>\r\n<\/pre>\n<\/div>\n
Collective communication<\/h2>\n
In parallel programming, it is often useful to perform what is known as collective communication, for example, broadcasting a Python object from the master process to all the worker processes. The example code below broadcasts a Numpy<\/a> array using the\u00a0bcast<\/code> method.<\/p>\n
\n
\n
from<\/span> mpi4py import<\/span> MPI\r\nimport<\/span> numpy as<\/span> np\r\n\r\ncomm =<\/span> MPI.<\/span>COMM_WORLD\r\nrank =<\/span> comm.<\/span>Get_rank(<\/span>)<\/span>\r\nsize =<\/span> comm.<\/span>Get_size(<\/span>)<\/span>\r\n\r\nif<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    data =<\/span> np.<\/span>arange(<\/span>4.0<\/span>)<\/span>\r\nelse<\/span>:<\/span>\r\n    data =<\/span> None<\/span>\r\n\r\ndata =<\/span> comm.<\/span>bcast(<\/span>data,<\/span> root=<\/span>0<\/span>)<\/span>\r\n\r\nif<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    print<\/span>(<\/span>'Process {} broadcast data:'<\/span>.<\/span>format(<\/span>rank)<\/span>,<\/span> data)<\/span>\r\nelse<\/span>:<\/span>\r\n    print<\/span>(<\/span>'Process {} received data:'<\/span>.<\/span>format(<\/span>rank)<\/span>,<\/span> data)<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n
The output from running the above code (broadcast.py) on 4 processes is as follows:<\/p>\n
\n
\n
$ srun -n 4 python3 broadcast.py\r\nProcess 0 broadcast data: [0. 1. 2. 3.]\r\nProcess 1 received data: [0. 1. 2. 3.]\r\nProcess 2 received data: [0. 1. 2. 3.]\r\nProcess 3 received data: [0. 1. 2. 3.]<\/pre>\n<\/div>\n<\/div>\n
If one needs to divide a task and distribute the sub-tasks to the processes, the scatter<\/code> method will be useful. However note that it is not possible to distribute more elements than the number of processors. If one has a big list or array, it is necessary to make slices of the list or array before calling the scatter<\/code> method. The code below is an example of scattering a Numpy<\/a> array, which is converted into a list of array slices before the scatter<\/code> method is called.<\/p>\n
\n
\n
from<\/span> mpi4py import<\/span> MPI\r\nimport<\/span> numpy as<\/span> np\r\n\r\ncomm =<\/span> MPI.<\/span>COMM_WORLD\r\nrank =<\/span> comm.<\/span>Get_rank(<\/span>)<\/span>\r\nnprocs =<\/span> comm.<\/span>Get_size(<\/span>)<\/span>\r\n\r\nif<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    data =<\/span> np.<\/span>arange(<\/span>15.0<\/span>)\r\n<\/span>\r\n    # determine the size of each sub-task<\/span>\r\n    ave,<\/span> res =<\/span> divmod<\/span>(<\/span>data.<\/span>size,<\/span> nprocs)<\/span>\r\n    counts =<\/span> [<\/span>ave +<\/span> 1<\/span> if<\/span> p <<\/span> res else<\/span> ave for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]\r\n<\/span>\r\n    # determine the starting and ending indices of each sub-task\r\n<\/span>    starts =<\/span> [<\/span>sum<\/span>(<\/span>counts[<\/span>:<\/span>p]<\/span>)<\/span> for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]<\/span>\r\n    ends =<\/span> [<\/span>sum<\/span>(<\/span>counts[<\/span>:<\/span>p+<\/span>1<\/span>]<\/span>)<\/span> for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]\r\n\r\n    # converts data into a list of arrays<\/span> <\/span>\r\n    data =<\/span> [<\/span>data[<\/span>starts[<\/span>p]<\/span>:<\/span>ends[<\/span>p]<\/span>]<\/span> for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]<\/span>\r\nelse<\/span>:<\/span>\r\n    data =<\/span> None<\/span>\r\n\r\ndata =<\/span> comm.<\/span>scatter(<\/span>data,<\/span> root=<\/span>0<\/span>)<\/span>\r\n\r\nprint<\/span>(<\/span>'Process {} has data:'<\/span>.<\/span>format(<\/span>rank)<\/span>,<\/span> data)<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n
The output from running the above code (scatter-array.py) on 4 processes is as follows.<\/p>\n
\n
\n
$ srun -n 4 python3 scatter-array.py\r\nProcess 0 has data: [0. 1. 2. 3.]\r\nProcess 1 has data: [4. 5. 6. 7.]\r\nProcess 2 has data: [8. 9. 10. 11.]\r\nProcess 3 has data: [12. 13. 14.]<\/pre>\n<\/div>\n<\/div>\n
Since there are 15 elements in total, we can see that each of the first 3 processes received 4 elements, and that the last process received 3 elements.<\/p>\n
The gather<\/code> method does the opposite to scatter<\/code>. If each process has an element, one can use gather<\/code> to collect them into a list of elements on the master process. The example code below uses scatter<\/code> and gather<\/code> to compute \u03c0 in parallel.<\/p>\n
\n
\n
from<\/span> mpi4py import<\/span> MPI\r\nimport<\/span> time\r\nimport<\/span> math\r\n\r\nt0 =<\/span> time.<\/span>time(<\/span>)<\/span>\r\n\r\ncomm =<\/span> MPI.<\/span>COMM_WORLD\r\nrank =<\/span> comm.<\/span>Get_rank(<\/span>)<\/span>\r\nnprocs =<\/span> comm.<\/span>Get_size(<\/span>)<\/span>\r\n\r\n# number of integration steps<\/span>\r\nnsteps =<\/span> 10000000<\/span>\r\n# step size<\/span>\r\ndx =<\/span> 1.0<\/span> \/<\/span> nsteps\r\n\r\nif<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    # determine the size of each sub-task<\/span>\r\n    ave,<\/span> res =<\/span> divmod<\/span>(<\/span>nsteps,<\/span> nprocs)<\/span>\r\n    counts =<\/span> [<\/span>ave +<\/span> 1<\/span> if<\/span> p <<\/span> res else<\/span> ave for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]\r\n<\/span>\r\n    # determine the starting and ending indices of each sub-task\r\n<\/span>    starts =<\/span> [<\/span>sum<\/span>(<\/span>counts[<\/span>:<\/span>p]<\/span>)<\/span> for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]<\/span>\r\n    ends =<\/span> [<\/span>sum<\/span>(<\/span>counts[<\/span>:<\/span>p+<\/span>1<\/span>]<\/span>)<\/span> for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]\r\n\r\n    # save the starting and ending indices in data <\/span> <\/span>\r\n    data =<\/span> [<\/span>(<\/span>starts[<\/span>p]<\/span>,<\/span> ends[<\/span>p]<\/span>)<\/span> for<\/span> p in<\/span> range<\/span>(<\/span>nprocs)<\/span>]<\/span>\r\nelse<\/span>:<\/span>\r\n    data =<\/span> None<\/span>\r\n\r\ndata =<\/span> comm.<\/span>scatter(<\/span>data,<\/span> root=<\/span>0<\/span>)<\/span>\r\n\r\n# compute partial contribution to pi on each process<\/span>\r\npartial_pi =<\/span> 0.0<\/span>\r\nfor<\/span> i in<\/span> range<\/span>(<\/span>data[<\/span>0<\/span>]<\/span>,<\/span> data[<\/span>1<\/span>]<\/span>)<\/span>:<\/span>\r\n    x =<\/span> (<\/span>i +<\/span> 0.5<\/span>)<\/span> *<\/span> dx\r\n    partial_pi +<\/span>=<\/span> 4.0<\/span> \/<\/span> (<\/span>1.0<\/span> +<\/span> x *<\/span> x)<\/span>\r\npartial_pi *<\/span>=<\/span> dx\r\n\r\npartial_pi =<\/span> comm.<\/span>gather(<\/span>partial_pi,<\/span> root=<\/span>0<\/span>)<\/span>\r\n\r\nif<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    print<\/span>(<\/span>'pi computed in {:.3f} sec'<\/span>.<\/span>format(<\/span>time.<\/span>time(<\/span>)<\/span> -<\/span> t0)<\/span>)<\/span>\r\n    print<\/span>(<\/span>'error is {}'<\/span>.<\/span>format(<\/span>abs<\/span>(<\/span>sum<\/span>(<\/span>partial_pi)<\/span> -<\/span> math.<\/span>pi)<\/span>)<\/span>)<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n
In addition to the gather<\/code> method which collects the elements into a list, one can also use the reduce<\/code> method to collect the results. The last five lines of the above example can be rewritten as follows, where \u201cop=MPI.SUM<\/code>\u201d requires that pi<\/code> is obtained as the sum of partial_pi<\/code> from all processes.<\/p>\n
\n
\n
pi =<\/span> comm.<\/span>reduce<\/span>(<\/span>partial_pi,<\/span> op=<\/span>MPI.<\/span>SUM<\/span>,<\/span> root=<\/span>0<\/span>)<\/span>\r\n\r\nif<\/span> rank ==<\/span> 0<\/span>:<\/span>\r\n    print<\/span>(<\/span>'pi computed in {:.3f} sec'<\/span>.<\/span>format(<\/span>time.<\/span>time(<\/span>)<\/span> -<\/span> t0)<\/span>)<\/span>\r\n    print<\/span>(<\/span>'error is {}'<\/span>.<\/span>format(<\/span>abs<\/span>(<\/span>pi -<\/span> math.<\/span>pi)<\/span>)<\/span>)<\/span>\r\n<\/pre>\n<\/div>\n<\/div>\n
Summary<\/h2>\n
We have briefly introduced the mpi4py<\/code> module and the ways it enables communication of Python objects via all-lowercase methods.<\/p>\n