A Random Walk with MPI
We would like to start with a simple example that has a simple work distribution. To make it even easier, we will use no communications.
Two-Dimensional Lattice Random Walk
A particle is moving through a two-dimensional finite grid. At each step, the “walker” can move left, right, up, or down, with equal probability. We seek the distance from the origin after N steps.
Theoretically, for N steps the distance is $\sqrt{N}$. We want to test this empirically, but one trial does not give us very good results. We want to run a lot of trials and average the results.
Agglomeration and Mapping
This algorithm has the properties that:
- Each trial is independent of the others
- There is a fixed number of tasks
- No communications is needed between tasks to obtain independent results.
Serial Code
Download the serial code for your language of choice. Compile it (if appropriate), then run a test case.
C++
/* Solves a random\-walk problem by a brute\-force method.
*
*/
#include <iostream>
#include <string>
#include <sstream>
#include <random>
#include <cmath>
using namespace std;
int main(int argc, char **argv) {
random_device rd;
mt19937 rng(rd());
uniform_int_distribution<int> choice(1,4);
long long N;
if (argc != 2) {
cout<<"0,0,0\n";
return 1;
}
else {
string steps=argv[1];
stringstream ssteps;
ssteps<<steps;
ssteps>>N;
}
double x=0;
double y=0;
int direction;
for (long long i=1; i<=N; ++i) {
direction=choice(rng);
switch (direction) {
case 1:
x+=1.;
break;
case 2:
x-=1.;
break;
case 3:
y+=1.;
break;
case 4:
y-=1.;
break;
}
}
double eucDist=sqrt(x*x+y*y);
cout<<N<<","<<sqrt((double)N)<<","<<eucDist<<"\n";
return 0;
}
Fortran
module random
implicit none
! Author: Katherine Holcomb
! Shortened version of a longer module
! Module incorporated into this file for convenience; typically it would be in its
! own file.
! Comment out one or the other when the module is incorporated into a code.
! Single precision
!integer, parameter :: rk = kind(1.0)
! Double precision
integer, parameter :: rk = kind(1.0d0)
contains
subroutine set_random_seed(seed)
! Sets all elements of the seed array
integer, optional, intent(in) :: seed
integer :: isize
integer, dimension(:), allocatable :: iseed
integer, dimension(8) :: idate
integer :: icount, i
call random_seed(size=isize)
if ( .not. present(seed) ) then
call system_clock(icount)
allocate(iseed(isize),source=icount+370*[(i,i=0,isize-1)])
else
allocate(iseed(isize))
iseed=seed
endif
call random_seed(put=iseed)
end subroutine set_random_seed
function urand(lb,ub,seed)
! Returns a uniformly-distributed random number in the range lb to ub.
real(rk) :: urand
real(rk), optional, intent(in) :: lb,ub
real(rk), optional, intent(in) :: seed
integer :: iseed
real(rk) :: rnd
real(rk) :: lower,upper
if ( present(seed) ) then
iseed=int(seed)
call set_random_seed(iseed)
endif
if ( present(lb) ) then
lower=lb
else
lower=0.0_rk
endif
if ( present(ub) ) then
upper = ub
else
upper = 1.0_rk
endif
call random_number(rnd)
urand = lower+(upper-lower)*rnd
return
end function urand
function randint(n,m)
! Returns a random integer between n and m.
integer :: randint
integer, intent(in) :: n,m
randint=ceiling(urand(real(n-1,rk),real(m,rk)))
end function randint
end module random
program random_walk
use random
implicit none
integer, parameter :: ik = selected_int_kind(15)
integer :: nargs
character(len=10) :: ns,nexp
integer :: m, step
integer(ik) :: n, nsteps
real, dimension(2) :: pos, move
real :: distance, avg
nargs=command_argument_count()
if ( nargs .ne. 1 ) then
write(*,*) 0, 0., 0.
stop
else
call get_command_argument(1,ns)
read(ns,'(i10)') nsteps
endif
call set_random_seed()
pos=[0.,0.]
avg=0.
do n=1,nsteps
step=randint(1,4)
if (step==1) then
move=[0.,1.]
else if (step==2) then
move=[0.,-1.]
else if (step==3) then
move=[1.,0.]
else
move=[-1.,0.]
endif
pos=pos+move
enddo
distance=sqrt(pos(1)**2+pos(2)**2)
write(*,*) nsteps, sqrt(real(nsteps)), distance
end program
Python
import numpy as np
import random
import sys
if len(sys.argv)>1:
Nsteps=int(sys.argv[1])
else:
print(0, 0., 0.)
sys.exit()
pos=np.array([0,0])
moves=[np.array([0,1]),np.array([0,-1]),np.array([1,0]),np.array([-1,0])]
for n in range(Nsteps):
move=random.choice(moves)
pos+=move
distance=np.sqrt(pos[0]**2+pos[1]**2)
print("{:.2f},{:.2f},{:.2f}".format(Nsteps, np.sqrt(Nsteps), distance))
Using the mpi1 program in your language as an example, add the lines to run this code with MPI. In your print statement make the first output value be the rank.
Run it
mpicxx -o randc mpirandom_walk.cxx
mpiexec -np 4 ./randc 1000000
0:1000000,1000,898.186
1:1000000,1000,189.589
2:1000000,1000,1235.87
3:1000000,1000,391.479
We have to compute the average manually
678.781
Try with -np 8
0:1000000,1000,1011.24
5:1000000,1000,1569.2
3:1000000,1000,1753.23
2:1000000,1000,967.042
1:1000000,1000,1212.17
6:1000000,1000,418.708
7:1000000,1000,881.862
4:1000000,1000,744.641
Why is the rank order jumbled?
MPI output is nondeterministic unless the programmer forces it to be ordered, using a barrier.