MPI Two-Dimensional Exchange
We can now use types to illustrate a two-dimensional halo exchange. For these examples, we will revert to blcoking send and receive to simplify the codes.
Special Note for Python
As we have mentioned, the internal values of a NumPy array are not directly accessibly by MPI. We can use the Vector type for column 0, but anything else rquires a different approach. One of the simpler methods is to use another built-in MPI type, the subarray. The basic syntax is
subarr=MPI.DOUBLE.Create_subarray(sizes,subsizes,starts,order=ORDER_C)
where sizes is an ordered type (such as a list) of the size along each dimension, subsizes is the size of the subarray to be selected, starts is the array of starting locations (zero based) for the subarray, and order is whether to use row-major (C) order or column-major (Fortran) order; ORDER_C is the default and is optional. In particular we can define
column_zero=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,0])
column_zero.Commit()
column_one=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,1])
column_one.Commit()
column_end=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,ncl])
column_end.Commit()
column_endbc=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,ncl+1])
column_endbc.Commit()
We then send and receive with syntax such as
comm.Recv([u,MPI.DOUBLE],source=MPI.ANY_SOURCE,tag=MPI.ANY_TAG,status=status)
comm.Send([w,MPI.DOUBLE],dest=0,tag=rank)
The subarray can also be used with C++ and Fortran, as we will describe in a later section.
Examples
Fortran
Contents of mpi_twod_exchange.cxx
#include <iostream>
#include <iomanip>
#include <cstring>
#include <mpi.h>
using namespace std;
/*#This example exchanges data among four rectangular domains with halos. Most real codes use squares, but we want to illustrate how to use different dimensions. Divide up the processes. Either we require a perfect square, or we must specify how to distribute by row/column. In a realistic program,
the process distribution (either the total, for a perfect square, or
the rows/columns) would be read in and we would need to check that the number
of processes requested is consistent with the decomposition. */
int main (int argc, char *argv[]) {
int N, M;
// Added for MPI
int nrl, ncl;
int rank, nprocs;
int up, down, right, left;
int tag=0;
//Initialize MPI
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
N = 8;
M = 12;
int nproc_rows=2;
int nproc_cols=3;
if (nproc_rows*nproc_cols != nprocs) {
cout<<"Number of rows times columns does not equal nprocs\n";
return 1;
}
if (N%nproc_rows==0 && M%nproc_cols==0) {
nrl = N/nproc_rows;
ncl = M/nproc_cols;
}
else {
cout<<"Number of ranks should divide the number of rows evenly.\n";
return 2;
}
double **w=new double*[nrl+2];
double *wptr=new double[(nrl+2)*(ncl+2)];
for (int i=0;i<nrl+2;++i,wptr+=ncl+2) {
w[i] = wptr;
}
double counter=1.;
for ( int i = 0; i < nrl+2; i++ ) {
for (int j = 0; j < ncl+2; j++ ) {
w[i][j] = (rank+1)*counter;
counter++;
}
}
//Set up the topology assuming processes number left to right by row
int my_row=rank/nproc_cols;
int my_col=rank%nproc_cols;
//Boundary conditions
double topBC=0.;
double bottomBC=200.;
double rightBC=100.;
double leftBC=100.;
if (my_col==0) {
for (int i=0;i<=nrl+1;++i){
w[i][0]=leftBC;
}
}
if (my_col==nproc_cols-1) {
for (int i=0;i<=nrl+1;++i){
w[i][ncl+1]=rightBC;
}
}
if (my_row==0) {
for (int i=0;i<=ncl+1;++i){
w[0][i]=topBC;
}
}
if (my_row==nproc_rows-1) {
for (int i=0;i<=ncl+1;++i){
w[nrl+1][i]=bottomBC;
}
}
//Find my neighbors
if (my_row==0) {
up=MPI_PROC_NULL;
}
else {
up=rank-nproc_cols;
}
if (my_row==nproc_rows-1) {
down=MPI_PROC_NULL;
}
else {
down= rank+nproc_cols;
}
if (my_col==0) {
left=MPI_PROC_NULL;
}
else {
left= rank-1;
}
if (my_col==nproc_cols-1) {
right=MPI_PROC_NULL;
}
else {
right = rank+1;
}
cout<<"Topology "<<rank<<" "<<up<<" "<<down<<" "<<left<<" "<<right<<endl;
MPI_Status status;
cout<<"Topology "<<rank<<" "<<up<<" "<<down<<" "<<left<<" "<<right<<endl;
//Force showing all of one rank array at a time. Not efficient.
int uwsize=(nrl+2)*(ncl+2);
double *u = new double[uwsize];
memset(u,0.,uwsize);
int position;
MPI_Barrier(MPI_COMM_WORLD);
if ( rank == 0 ) {
cout<<"---W Before for rank 0---"<<endl;
for (int i=0;i<nrl+2;i++) {
for (int j=0;j<ncl+2;j++) {
cout<<w[i][j]<<" ";
}
cout<<endl;
}
for (int n=1;n<nprocs;++n) {
memset(u,0.,uwsize);
MPI_Recv(u,uwsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
cout<<"--Before for rank "<<status.MPI_SOURCE<<endl;
position=0;
for (int i=0;i<=nrl+1;i++) {
for (int j=0;j<=ncl+1;j++) {
cout<<u[position++]<<" ";
}
cout<<endl;
}
}
}
else {
// Pack the 2D array into the buffer
position=0;
for (int i=0; i<=nrl+1; i++)
for (int j=0;j<=ncl+1; j++)
u[position++]=w[i][j];
MPI_Send(u,uwsize, MPI_DOUBLE, 0, rank, MPI_COMM_WORLD);
}
MPI_Barrier(MPI_COMM_WORLD);
//Set up the MPI type for columns
MPI_Datatype col;
MPI_Type_vector(nrl,1,ncl+2,MPI_DOUBLE,&col);
MPI_Type_commit(&col);
MPI_Sendrecv(&w[1][0],ncl+2, MPI_DOUBLE,up,tag,&w[nrl+1][0],
ncl+2, MPI_DOUBLE,down,tag,MPI_COMM_WORLD,&status);
MPI_Sendrecv(&w[nrl][0],ncl+2,MPI_DOUBLE,down,tag,&w[0][0],
ncl+2,MPI_DOUBLE,up,tag,MPI_COMM_WORLD,&status);
MPI_Sendrecv(&w[1][ncl],1,col,right,tag,&w[1][0],1,col,left, tag, MPI_COMM_WORLD, &status);
MPI_Sendrecv(&w[1][1],1,col,left,tag,&w[1][ncl+1],1,col,right, tag, MPI_COMM_WORLD, &status);
//Print results
MPI_Barrier(MPI_COMM_WORLD);
if ( rank == 0 ) {
cout<<"---After for rank 0---"<<endl;
for (int i=0;i<nrl+2;i++) {
for (int j=0;j<ncl+2;j++) {
cout<<w[i][j]<<" ";
}
cout<<endl;
}
for (int n=1;n<nprocs;++n) {
memset(u,0.,uwsize);
MPI_Recv(u,uwsize,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
cout<<"--After for rank "<<status.MPI_SOURCE<<endl;
position=0;
for (int i=0;i<=nrl+1;i++) {
for (int j=0;j<=ncl+1;j++) {
cout<<u[position++]<<" ";
}
cout<<endl;
}
}
}
else {
// Pack the 2D array into the buffer
position=0;
for (int i=0; i<=nrl+1; i++)
for (int j=0;j<=ncl+1; j++)
u[position++]=w[i][j];
MPI_Send(u,uwsize, MPI_DOUBLE, 0, rank, MPI_COMM_WORLD);
}
MPI_Barrier(MPI_COMM_WORLD);
//end checking section
//Type_free not really necessary but good practice
MPI_Type_free(&col);
MPI_Finalize();
return 0;
}
Fortran
Contents of mpi_twod_exchange.f90
program twod_exchange
use mpi_f08
implicit none
integer :: N, M
integer :: i,j
integer :: rank, nprocs, tag=0
integer :: nproc_rows, nproc_cols, my_row, my_col, nrl, ncl
integer :: left, right, up, down
integer :: uwsize
double precision, dimension(:,:), allocatable :: w, u
double precision :: topBC, bottomBC, leftBC, rightBC
type(MPI_Status) :: mpi_stat
type(MPI_Datatype) :: row
!Initialize MPI, get the local number of columns
call MPI_INIT()
call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs)
call MPI_COMM_RANK(MPI_COMM_WORLD,rank)
!This example exchanges data among four rectangular domains with halos.
!Most real codes use squares, but we want to illustrate how to use different
!dimensions.
!Divide up the processes. Either we require a perfect square, or we
!must specify how to distribute by row/column. In a realistic program,
!the process distribution (either the total, for a perfect square, or
!the rows/columns) would be read in and we would need to check that the number
!of processes requested is consistent with the decomposition.
N = 8
M = 12
nproc_rows=2
nproc_cols=3
if (nproc_rows*nproc_cols /= nprocs) then
call MPI_Finalize()
stop "Number of rows times columns does not equal nprocs"
endif
if (mod(N,nproc_rows)==0 .and. mod(M,nproc_cols)==0) then
nrl = N/nproc_rows
ncl = M/nproc_cols
else
call MPI_Finalize()
stop "Number of ranks should divide the number of rows evenly."
endif
allocate(w(0:nrl+1, 0:ncl+1),u(0:nrl+1,0:ncl+1))
w=reshape([(i,i=1,(nrl+2)*(ncl+2))],(/nrl+2,ncl+2/))
w=(rank+1)*w
!Set up the topology assuming processes numbered left to right by row
my_row=rank/nproc_cols
my_col=mod(rank,nproc_cols)
! Setting up the neighbors for each process
if (my_row == 0) then
up = MPI_PROC_NULL
else
up = rank - nproc_cols
endif
if (my_row == nproc_rows-1) then
down = MPI_PROC_NULL
else
down = rank + nproc_cols
endif
if (my_col == 0) then
left = MPI_PROC_NULL
else
left = rank-1
endif
if (my_col == nproc_cols-1) then
right = MPI_PROC_NULL
else
right = rank+1
endif
! Boundary conditions
topBC=0.d0
bottomBC=200.d0
leftBC=300.d0
rightBC=300.d0
if (my_row==0) then
w(0,:)=topBC
endif
if (my_row==nproc_rows-1) then
w(nrl+1,:)=bottomBC
endif
if (my_col==0) then
w(:,0)=leftBC
endif
if (my_col==nproc_cols-1) then
w(:,ncl+1)=rightBC
endif
call MPI_Type_vector(ncl+2, 1, nrl+2, MPI_DOUBLE_PRECISION, row)
call MPI_Type_commit(row)
tag=0
write(*,*) "Topology ",rank, left, right
!starting
!This forces the output to show one rank at a time. It's not efficient.
uwsize=size(u)
if (rank==0) then
write(*,*) 'Initial for rank 0'
do j=0,nrl+1
write(*,'(*(f10.2))') w(j,:)
enddo
do i=1,nprocs-1
u=0.d0
call MPI_Recv(u,uwsize,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE, &
MPI_ANY_TAG, MPI_COMM_WORLD,mpi_stat)
write(*,*) 'Initial for rank',mpi_stat%MPI_SOURCE
do j=0,nrl+1
write(*,'(*(f10.2))') u(j,:)
enddo
enddo
else
call MPI_Send(w,uwsize,MPI_DOUBLE_PRECISION,0,rank,MPI_COMM_WORLD)
endif
call MPI_Barrier(MPI_COMM_WORLD)
! Exchange halo values
! Send left and right
call MPI_SENDRECV(w(1:nrl,1) ,nrl,MPI_DOUBLE_PRECISION,left,tag, &
w(1:nrl,ncl+1),nrl,MPI_DOUBLE_PRECISION,right,tag, &
MPI_COMM_WORLD,mpi_stat)
call MPI_SENDRECV(w(1:nrl,ncl) ,nrl,MPI_DOUBLE_PRECISION,right,tag, &
w(1:nrl,0) ,nrl,MPI_DOUBLE_PRECISION,left,tag, &
MPI_COMM_WORLD,mpi_stat)
! Send up and down
call MPI_SENDRECV(w(1,0),1, row, up,tag, w(nrl+1,0),1, row,down,tag, &
MPI_COMM_WORLD,mpi_stat)
call MPI_SENDRECV(w(nrl,0),1,row, down,tag, w(0,0) ,1,row, up,tag, &
MPI_COMM_WORLD,mpi_stat)
!result
if (rank==0) then
write(*,*) 'Result for rank 0'
do j=0,nrl+1
write(*,'(*(f10.2))') w(j,:)
enddo
do i=1,nprocs-1
u=0.d0
call MPI_Recv(u,uwsize,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE, &
MPI_ANY_TAG, MPI_COMM_WORLD,mpi_stat)
write(*,*) 'Result for rank',mpi_stat%MPI_SOURCE
do j=0,nrl+1
write(*,'(*(f10.2))') u(j,:)
enddo
enddo
else
call MPI_Send(w,uwsize,MPI_DOUBLE_PRECISION,0,rank,MPI_COMM_WORLD)
endif
call MPI_Type_free(row)
call MPI_Finalize()
end program
Python
Contents of mpi_twod_exchange.py
import sys
import numpy as np
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
nprocs = comm.Get_size()
#This example exchanges data among four rectangular domains with halos.
#Most real codes use squares, but we want to illustrate how to use different
#dimensions.
#Divide up the processes. Either we require a perfect square, or we
#must specify how to distribute by row/column. In a realistic program,
#the process distribution (either the total, for a perfect square, or
#the rows/columns) would be read in and we would need to check that the number
#of processes requested is consistent with the decomposition.
N = 8
M = 12
nproc_rows=2
nproc_cols=3
if nproc_rows*nproc_cols != nprocs:
print("Number of rows times columns does not equal nprocs")
sys.exit()
if N%nproc_rows==0 and M%nproc_cols==0:
nrl = N//nproc_rows
ncl = M//nproc_cols
else:
print("Number of ranks should divide the number of rows evenly.")
sys.exit()
w = np.zeros((nrl+2, ncl+2), dtype=np.double)
#Set up the topology assuming processes numbered left to right by row
my_row=rank//nproc_cols
my_col=rank%nproc_cols
#Arbitrary values.
w[:,:] = np.reshape(np.arange(1,(nrl+2)*(ncl+2)+1),(nrl+2,ncl+2))
w=w*(rank+1)
# setting up the up and down rank for each process
if my_row == 0 :
up = MPI.PROC_NULL
else :
up = rank - nproc_cols
if my_row == nproc_rows-1 :
down = MPI.PROC_NULL
else :
down = rank + nproc_cols
if my_col == 0 :
left = MPI.PROC_NULL
else:
left = rank-1
if my_col == nproc_cols-1:
right = MPI.PROC_NULL
else:
right = rank+1
#Boundary conditions
topBC=0.0
bottomBC=200.00
rightBC=100.0
leftBC=100.0
if right == MPI.PROC_NULL:
w[:,ncl+1]=rightBC
if left == MPI.PROC_NULL:
w[:,0]=leftBC
if up == MPI.PROC_NULL:
w[0,:]=topBC
if down == MPI.PROC_NULL:
w[nrl+1,:]=bottomBC
# set up MPI type for left column
column_zero=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,0])
column_zero.Commit()
column_one=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,1])
column_one.Commit()
column_end=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,ncl])
column_end.Commit()
column_endbc=MPI.DOUBLE.Create_subarray([nrl+2,ncl+2],[nrl,1],[1,ncl+1])
column_endbc.Commit()
#This forces the output to show one one rank at a time. It is not efficient.
status=MPI.Status()
print("topology", rank, up, down, left, right)
if rank==0:
u=np.zeros_like(w)
print("Initial for rank 0")
print(w)
for i in range(1,nprocs):
comm.Recv([u,MPI.DOUBLE],source=MPI.ANY_SOURCE,tag=MPI.ANY_TAG,status=status)
print("Initial for rank",status.Get_source())
print(u)
else:
comm.Send([w,MPI.DOUBLE],dest=0,tag=rank)
comm.Barrier()
tag=0
# sending up and receiving down
comm.Sendrecv([w[1,0:ncl+2],MPI.DOUBLE], up, tag, [w[nrl+1,0:ncl+2],MPI.DOUBLE], down, tag)
# sending down and receiving up
comm.Sendrecv([w[nrl,0:ncl+2],MPI.DOUBLE], down, tag, [w[0,0:ncl+2],MPI.DOUBLE], up, tag)
# sending left and right
comm.Sendrecv([w,1,column_one], left, tag, [w,1,column_endbc], right, tag)
comm.Sendrecv([w,1,column_end], right, tag, [w,1,column_zero], left, tag)
# check result
if rank==0:
u=np.zeros_like(w)
print("Final for rank 0")
print(w)
for i in range(1,nprocs):
comm.Recv([u,MPI.DOUBLE],source=MPI.ANY_SOURCE,tag=MPI.ANY_TAG,status=status)
print("Final for rank",status.Get_source())
print(u)
else:
comm.Send([w,MPI.DOUBLE],dest=0,tag=rank)
# Wrap up
column_zero.Free()
column_one.Free()
column_end.Free()
column_endbc.Free()