Introduction to Parallel Programming

Probability

${\displaystyle f:U\subset \mathbb {R} ^{2}\rightarrow \mathbb {R} _{+}}$ is a probability density function if ${\displaystyle \int \int _{U}f\mathrm {d} A=1}$

If ${\displaystyle f}$ is a probability density function which takes the set ${\displaystyle U\subset \mathbb {R} ^{2}}$, then the probability of events in the set ${\displaystyle W\subset U}$ occurring is ${\displaystyle P(W)=\int \int _{W}f\mathrm {d} A.}$

The joint density for it to snow ${\displaystyle x}$ inches tomorrow and for Kelly to win ${\displaystyle y}$ dollar in the lottery tomorrow is given by ${\displaystyle f={\frac {c}{(1+x)(100+y)}}}$ for ${\displaystyle x,y\in [0,100]\times [0,100]}$ and ${\displaystyle f=0}$ otherwise. Find ${\displaystyle c}$.

Suppose ${\displaystyle X}$ is a random variable with probability density function ${\displaystyle f_{1}(x)}$ and ${\displaystyle Y}$ is a random variable with a probability density function ${\displaystyle f_{2}(y)}$. Then ${\displaystyle X}$ and ${\displaystyle Y}$ are independent random variables if their joint density function is ${\displaystyle f(x,y)=f_{1}(x)f_{2}(y).}$

The probability it will snow tomorrow and the probability Kelly will win the lottery tomorrow are independent random variables.

If ${\displaystyle f(x,y)}$ is a probability density function for the random variables ${\displaystyle X}$ and ${\displaystyle Y}$, the X mean is ${\displaystyle \mu _{1}={\bar {X}}=\int \int xf\mathrm {d} A}$ and the Y mean is ${\displaystyle \mu _{2}={\bar {Y}}=\int \int yf\mathrm {d} A.}$

The X mean and the Y mean are the expected values of X and Y.

If ${\displaystyle f(x,y)}$ is a probability density function for the random variables ${\displaystyle X}$ and ${\displaystyle Y}$, the X variance is ${\displaystyle \sigma _{1}^{2}={\overline {(X-{\bar {X}})^{2}}}=\int \int (x-{\bar {X}})^{2}f\mathrm {d} A}$ and the Y variance is ${\displaystyle \sigma _{2}^{2}={\overline {(Y-{\bar {Y}})^{2}}}=\int \int (y-{\bar {Y}})^{2}f\mathrm {d} A.}$

The standard deviation is defined to be the square root of the variance.

Find an expression for the probability that it will snow more than 1.1 times the expected snowfall and also that Kelly will win more than 1.2 times the expected amount in the lottery.

Exercise

• A class is graded on a curve. It is assumed that the class is a representative sample of the population, the probability density function for the numerical score ${\displaystyle x}$ is given by ${\displaystyle f(x)=C\exp \left(-{\frac {(x-\mu )^{2}}{2\sigma ^{2}}}\right).}$ For simplicity we assume that ${\displaystyle x}$ can take on the values ${\displaystyle -\infty }$ and ${\displaystyle \infty }$, though in actual fact the exam is scored from ${\displaystyle 0}$ to ${\displaystyle 100}$.

  1. Determine ${\displaystyle C}$ using results from your previous homework.

  2. Suppose there are 240 students in the class and the mean and standard deviation for the class is not reported. As an enterprising student, you poll 60 of your fellow students (we shall suppose they are selected randomly). You find that the mean for these 60 students is 55% and the standard deviation is 10%. Use the Student’s t distribution http://en.wikipedia.org/wiki/Student\%27s_t-distribution to estimate the 90% confidence interval for the actual sample mean. Make a sketch of the t-distribution probability density function and shade the region which corresponds to the 90% confidence interval for the sample mean.^[5]

  Remark Fortunately, all the students are hard working, so the possibility of a negative score, although possible, is extremely low, and so we neglect it to make the above computation easier.

Riemann Integration

Recall that we can approximate integrals by Riemann sums. There are many integrals one cannot evaluate analytically, but for which a numerical answer is required. In this section, we shall explore a simple way of doing this on a computer. Suppose we want to find ${\displaystyle I2d=\int _{0}^{1}\int _{0}^{4}x^{2}+2y^{2}\mathrm {d} y\mathrm {d} x.}$ If we do this analytically we find ${\displaystyle I2d=44.}$ Let us suppose we have forgotten how to integrate, and so we do this numerically. We can do so using the following Matlab code:

% A program to approximate an integral

clear all; format compact; format short;

nx=1000; % number of points in x
xend=1; % last discretization point
xstart=0; % first discretization point
dx=(xend-xstart)/(nx-1); % size of each x sub-interval

ny=4000; % number of points in y
yend=4; % last discretization point
ystart=0; % first discretization point
dy=(yend-ystart)/(ny-1); % size of each y sub-interval

% create vectors with points for x and y
for i=1:nx
    x(i)=xstart+(i-1)*dx;
end
for j=1:ny
    y(j)=ystart+(j-1)*dy;
end

% Approximate the integral by a sum
I2d=0;
for i=1:nx
    for j=1:ny
        I2d=I2d+(x(i)^2+2*y(j)^2)*dy*dx;
    end
end
% print out final answer
I2d

We can do something similar in three dimensions. Suppose we want to calculate ${\displaystyle I3d=\int _{0}^{1}\int _{0}^{1}\int _{0}^{4}x^{2}+2y^{2}+3z^{2}\mathrm {d} z\mathrm {d} y\mathrm {d} x.}$ Analytically we find that ${\displaystyle I3d=68}$

Exercises

1. Modify the Matlab code to perform the three dimensional integral.
2. Try and determine how the accuracy of either the two or three dimensional method varies as the number of subintervals is changed.

Monte Carlo Integration

It is possible to extend the above integration schemes to higher and higher dimensional integrals. This can become computationally intensive and an alternate method of integration based on probability is often used. The method we will discuss is called the Monte Carlo method, and the description of this method which follows is based on Chapter 3 of Vector Calculus by Michael Corral^[6] which is available at http://www.mecmath.net/ and where Java and Sage programs for doing Monte Carlo integration can be found. The idea behind Monte Carlo integration is based on the concept of the average value of a function, which is encountered in single-variable calculus. Recall that for a continuous function ${\displaystyle f(x)}$, the average value ${\displaystyle {\bar {f}}}$ of ${\displaystyle f}$ over an interval ${\displaystyle \lbrack a,b\rbrack }$ is defined as

${\displaystyle {\bar {f}}~=~{\frac {1}{b-a}}\int _{a}^{b}f(x)\,dx~.}$

( 1)

The quantity ${\displaystyle b-a}$ is the length of the interval ${\displaystyle \lbrack a,b\rbrack }$, which can be thought of as the “volume” of the interval. Applying the same reasoning to functions of two or three variables, we define the average value of ${\displaystyle f(x,y)}$ over a region ${\displaystyle R}$ to be

${\displaystyle {\bar {f}}~=~{\frac {1}{A(R)}}\iint \limits _{R}f(x,y)\,dA~,}$

( 2)

where ${\displaystyle A(R)}$ is the area of the region ${\displaystyle R}$, and we define the average value of ${\displaystyle f(x,y,z)}$ over a solid ${\displaystyle S}$ to be

${\displaystyle {\bar {f}}~=~{\frac {1}{V(S)}}\iiint \limits _{S}f(x,y,z)\,dV~,}$

( 3)

where ${\displaystyle V(S)}$ is the volume of the solid ${\displaystyle S}$. Thus, for example, we have

${\displaystyle \iint \limits _{R}f(x,y)\,dA~=~A(R){\bar {f}}~.}$

( 4)

The average value of ${\displaystyle f(x,y)}$ over ${\displaystyle R}$ can be thought of as representing the sum of all the values of ${\displaystyle f}$ divided by the number of points in ${\displaystyle R}$. Unfortunately there are an infinite number (in fact, uncountably many) points in any region, i.e. they can not be listed in a discrete sequence. But what if we took a very large number ${\displaystyle N}$ of random points in the region ${\displaystyle R}$ (which can be generated by a computer) and then took the average of the values of ${\displaystyle f}$ for those points, and used that average as the value of ${\displaystyle {\bar {f}}}$? This is exactly what the Monte Carlo method does. So in formula 4 the approximation we get is

${\displaystyle \iint \limits _{R}f(x,y)\,dA~\approx ~A(R){\bar {f}}\pm A(R){\sqrt {\frac {{\overline {f^{2}}}-({\bar {f}})^{2}}{N}}}~,}$

( 5)

where

${\displaystyle {\bar {f}}~=~{\frac {\sum _{i=1}^{N}f(x_{i},y_{i})}{N}}\quad {\text{and}}\quad {\overline {f^{2}}}~=~{\frac {\sum _{i=1}^{N}(f(x_{i},y_{i}))^{2}}{N}}~,}$

( 6)

with the sums taken over the ${\displaystyle N}$ random points ${\displaystyle (x_{1},y_{1})}$, ${\displaystyle \ldots }$, ${\displaystyle (x_{N},y_{N})}$. The ${\displaystyle \pm }$ “error term” in formula 5 does not really provide hard bounds on the approximation. It represents a single standard deviation from the expected value of the integral. That is, it provides a likely bound on the error. Due to its use of random points, the Monte Carlo method is an example of a probabilistic method (as opposed to deterministic methods such as the Riemann sum approximation method, which use a specific formula for generating points).

For example, we can use the formula in eq. 5 to approximate the volume ${\displaystyle V}$ under the surface ${\displaystyle z=x^{2}+2y^{2}}$ over the rectangle ${\displaystyle R=(0,1)\times (0,4)}$. Recall that the actual volume is ${\displaystyle 44}$. Below is a Matlab code that calculates the volume using Monte Carlo integration

( H)

A Matlab program which demonstrates how to use the Monte Carlo method to calculate the volume below ${\displaystyle z=x^{2}+2y^{2}}$, with ${\displaystyle (x,y)\in (0,1)\times (0,4)}$. Code download

% A program to approximate an integral using the Monte Carlos method

% This program can be made much faster by using Matlab's matrix and vector
% operations, however to allow easy translation to other languages we have
% made it as simple as possible.

Numpoints=65536; % number of random points

I2d=0; % Initialize value
I2dsquare=0; % initial variance
for n=1:Numpoints
    % generate random number drawn from a uniform distribution on (0,1)
    x=rand(1);
    y=rand(1)*4;
    I2d=I2d+x^2+2*y^2;
    I2dsquare=I2dsquare+(x^2+2*y^2)^2;
end
% we scale the integral by the total area and divide by the number of
% points used
I2d=I2d*4/Numpoints
% we also output an estimated error
I2dsquare=I2dsquare*4/Numpoints;
EstimError=4*sqrt( (I2d^2-I2dsquare)/Numpoints)

The results of running this program with various numbers of random points are shown below:

N = 16: 41.3026 +/- 30.9791
N = 256: 47.1855 +/- 9.0386
N = 4096: 43.4527 +/- 2.0280
N = 65536: 44.0026 +/- 0.5151

As you can see, the approximation is fairly good. As ${\displaystyle N\to \infty }$, it can be shown that the Monte Carlo approximation converges to the actual volume (on the order of ${\displaystyle O({\sqrt {N}})}$, in computational complexity terminology).

In the above example the region ${\displaystyle R}$ was a rectangle. To use the Monte Carlo method for a nonrectangular (bounded) region ${\displaystyle R}$, only a slight modification is needed. Pick a rectangle ${\displaystyle {\tilde {R}}}$ that encloses ${\displaystyle R}$, and generate random points in that rectangle as before. Then use those points in the calculation of ${\displaystyle {\bar {f}}}$ only if they are inside ${\displaystyle R}$. There is no need to calculate the area of ${\displaystyle R}$ for formula ({eqn:monte}) in this case, since the exclusion of points not inside ${\displaystyle R}$ allows you to use the area of the rectangle ${\displaystyle {\tilde {R}}}$ instead, similar to before.

For instance, one can show that the volume under the surface ${\displaystyle z=1}$ over the nonrectangular region ${\displaystyle R=\left\{(x,y):0\leq x^{2}+y^{2}\leq 1\right\}}$ is ${\displaystyle \pi }$. Since the rectangle ${\displaystyle {\tilde {R}}=[-1,1]\times [-1,1]}$ contains ${\displaystyle R}$, we can use a similar program to the one we used, the largest change being a check to see if ${\displaystyle y^{2}+x^{2}\leq 1}$ for a random point ${\displaystyle (x,y)}$ in ${\displaystyle [-1,1]\times [-1,1]}$. A Matlab code listing which demonstrates this is below:

( I)

A Matlab program which demonstrates how to use the Monte Carlo method to calculate the area of an irregular region and also to calculate ${\displaystyle \pi }$. Code download

% A program to approximate an integral using the Monte Carlos method

% This program can be made much faster by using Matlab's matrix and vector
% operations, however to allow easy translation to other languages we have
% made it as simple as possible.

Numpoints=256; % number of random points

I2d=0; % Initialize value
I2dsquare=0; % initial variance
for n=1:Numpoints
    % generate random number drawn from a uniform distribution on (0,1) and
    % scale this to (-1,1)
    x=2*rand(1)-1;
    y=2*rand(1) -1;
    if ((x^2+y^2) <1)
        I2d=I2d+1;
        I2dsquare=I2dsquare+1;
    end
end
% We scale the integral by the total area and divide by the number of
% points used
I2d=I2d*4/Numpoints
% we also output an estimated error
I2dsquare=I2dsquare*4/Numpoints;
EstimError=4*sqrt( (I2d^2-I2dsquare)/Numpoints)

The results of running the program with various numbers of random points are shown below:

N = 16: 3.5000 +/- 2.9580
N = 256: 3.2031 +/- 0.6641
N = 4096: 3.1689 +/- 0.1639
N = 65536: 3.1493 +/- 0.0407

To use the Monte Carlo method to evaluate triple integrals, you will need to generate random triples ${\displaystyle (x,y,z)}$ in a parallelepiped, instead of random pairs ${\displaystyle (x,y)}$ in a rectangle, and use the volume of the parallelepiped instead of the area of a rectangle in formula 5 . For a more detailed discussion of numerical integration methods, please take a further course in mathematics.

Exercises

1. Write a program that uses the Monte Carlo method to approximate the double integral ${\displaystyle \iint \limits _{R}e^{xy}\,dA}$, where ${\displaystyle R=\lbrack 0,1\rbrack \times \lbrack 0,1\rbrack }$. Show the program output for ${\displaystyle N=10}$, ${\displaystyle 100}$, ${\displaystyle 1000}$, ${\displaystyle 10000}$, ${\displaystyle 100000}$ and ${\displaystyle 1000000}$ random points.
2. Write a program that uses the Monte Carlo method to approximate the triple integral ${\displaystyle \iiint \limits _{S}e^{xyz}\,dV}$, where ${\displaystyle S=\lbrack 0,1\rbrack \times \lbrack 0,1\rbrack \times \lbrack 0,1\rbrack }$. Show the program output for ${\displaystyle N=10}$, ${\displaystyle 100}$, ${\displaystyle 1000}$, ${\displaystyle 10000}$, ${\displaystyle 100000}$ and ${\displaystyle 1000000}$ random points.
3. Use the Monte Carlo method to approximate the volume of a sphere of radius ${\displaystyle 1}$.

Parallel Monte Carlo Integration

As you may have noticed, the algorithms are simple, but can require very many grid points to become accurate. It is therefore useful to run these algorithms on a parallel computer. We will demonstrate a parallel Monte Carlo calculation of ${\displaystyle \pi }$. Before we can do this, we need to learn how to use a parallel computer^[7].

We now examine a Fortran program for calculating ${\displaystyle \pi }$. These programs are taken from http://chpc.wustl.edu/mpi-fortran.html, where further explanation can be found. The original source of these programs appears to be Gropp, Lusk and Skjellum^[8]

Serial

( I)

A serial Fortran program which demonstrates how to calculate ${\displaystyle \pi }$ using a Monte Carlo method Code download

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program use a monte carlo method to calculate pi
!
! .. Parameters ..
! npts = total number of Monte Carlo points
! xmin = lower bound for integration region
! xmax = upper bound for integration region
! .. Scalars ..
! i = loop counter
! f = average value from summation
! sum = total sum
! randnum = random number generated from (0,1) uniform
! distribution
! x = current Monte Carlo location
! .. Arrays ..
!
! .. Vectors ..
!
! REFERENCES
! http://chpc.wustl.edu/mpi-fortran.html
! Gropp, Lusk and Skjellum, "Using MPI" MIT press (1999)
!
! ACKNOWLEDGEMENTS
! The program below was modified from one available at the internet
! address in the references. This internet address was last checked
! on 30 March 2012
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
!
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! None
  PROGRAM monte_carlo
    IMPLICIT NONE

INTEGER(kind=8), PARAMETER :: npts = 1e10
    REAL(kind=8), PARAMETER :: xmin=0.0d0,xmax=1.0d0
    INTEGER(kind=8) :: i
    REAL(kind=8) :: f,sum, randnum,x

    DO i=1,npts
      CALL random_number(randnum)
      x = (xmax-xmin)*randnum + xmin
      sum = sum + 4.0d0/(1.0d0 + x**2)
    END DO
f = sum/npts
    PRINT*,'PI calculated with ',npts,' points = ',f

    STOP
END

( J)

An example makefile for compiling the program in listing I Code download

#define the complier
COMPILER = mpif90
# compilation settings, optimization, precision, parallelization
FLAGS = -O0

# libraries
LIBS =
# source list for main program
SOURCES = montecarloserial.f90

test: $(SOURCES)
${COMPILER} -o montecarloserial $(FLAGS) $(SOURCES)

clean:
rm *.o

clobber:
rm montecarloserial

( K)

An example submission script for use on Trestles located at the San Diego Supercomputing Center Code download

#!/bin/bash
# the queue to be used.
#PBS -q shared
# specify your project allocation
#PBS -A mia122
# number of nodes and number of processors per node requested
#PBS -l nodes=1:ppn=1
# requested Wall-clock time.
#PBS -l walltime=00:05:00
# name of the standard out file to be "output-file".
#PBS -o job_output
# name of the job
#PBS -N MCserial
# Email address to send a notification to, change "youremail" appropriately
#PBS -M youremail@umich.edu
# send a notification for job abort, begin and end
#PBS -m abe
#PBS -V
cd $PBS_O_WORKDIR #change to the working directory
mpirun_rsh -np 1 -hostfile $PBS_NODEFILE montecarloserial

Parallel

( L)

A parallel fortran MPI program to calculate ${\displaystyle \pi }$. Code download

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program uses MPI to do a parallel monte carlo calculation of pi
!
! .. Parameters ..
! npts = total number of Monte Carlo points
! xmin = lower bound for integration region
! xmax = upper bound for integration region
! .. Scalars ..
! mynpts = this processes number of Monte Carlo points
! myid = process id
! nprocs = total number of MPI processes
! ierr = error code
! i = loop counter
! f = average value from summation
! sum = total sum
! mysum = sum on this process
! randnum = random number generated from (0,1) uniform
! distribution
! x = current Monte Carlo location
! start = simulation start time
! finish = simulation end time
! .. Arrays ..
!
! .. Vectors ..
!
! REFERENCES
! http://chpc.wustl.edu/mpi-fortran.html
! Gropp, Lusk and Skjellum, "Using MPI" MIT press (1999)
!
! ACKNOWLEDGEMENTS
! The program below was modified from one available at the internet
! address in the references. This internet address was last checked
! on 30 March 2012
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
!
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! MPI library
    PROGRAM monte_carlo_mpi
    USE MPI
    IMPLICIT NONE

    INTEGER(kind=8), PARAMETER :: npts = 1e10
    REAL(kind=8), PARAMETER :: xmin=0.0d0,xmax=1.0d0
    INTEGER(kind=8) :: mynpts
    INTEGER(kind=4) :: ierr, myid, nprocs
    INTEGER(kind=8) :: i
    REAL(kind=8) :: f,sum,mysum,randnum
    REAL(kind=8) :: x, start, finish
    
    ! Initialize MPI
    CALL MPI_INIT(ierr)
    CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
    CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)
    start=MPI_WTIME()

! Calculate the number of points each MPI process needs to generate
    IF (myid .eq. 0) THEN
mynpts = npts - (nprocs-1)*(npts/nprocs)
    ELSE
mynpts = npts/nprocs
    ENDIF
    
    ! set initial sum to zero
    mysum = 0.0d0
! use loop on local process to generate portion of Monte Carlo integral
    DO i=1,mynpts
      CALL random_number(randnum)
      x = (xmax-xmin)*randnum + xmin
      mysum = mysum + 4.0d0/(1.0d0 + x**2)
    ENDDO

! Do a reduction and sum the results from all processes
    CALL MPI_REDUCE(mysum,sum,1,MPI_DOUBLE_PRECISION,MPI_SUM,&
          0,MPI_COMM_WORLD,ierr)
    finish=MPI_WTIME()

    ! Get one process to output the result and running time
    IF (myid .eq. 0) THEN
f = sum/npts
         PRINT*,'PI calculated with ',npts,' points = ',f
         PRINT*,'Program took ', finish-start, ' for Time stepping'
    ENDIF

CALL MPI_FINALIZE(ierr)

    STOP
END PROGRAM

( M)

An example makefile for compiling the program in listing L Code download

#define the complier
COMPILER = mpif90
# compilation settings, optimization, precision, parallelization
FLAGS = -O0

# libraries
LIBS =
# source list for main program
SOURCES = montecarloparallel.f90

test: $(SOURCES)
${COMPILER} -o montecarloparallel $(FLAGS) $(SOURCES)

clean:
rm *.o

clobber:
rm montecarloparallel

( N)

An example submission script for use on Trestles located at the San Diego Supercomputing Center Code download

#!/bin/bash
# the queue to be used.
#PBS -q normal
# specify your project allocation
#PBS -A mia122
# number of nodes and number of processors per node requested
#PBS -l nodes=1:ppn=32
# requested Wall-clock time.
#PBS -l walltime=00:05:00
# name of the standard out file to be "output-file".
#PBS -o job_output
# name of the job, you may want to change this so it is unique to you
#PBS -N MPI_MCPARALLEL
# Email address to send a notification to, change "youremail" appropriately
#PBS -M youremail@umich.edu
# send a notification for job abort, begin and end
#PBS -m abe
#PBS -V

# change to the job submission directory
cd $PBS_O_WORKDIR
# Run the job
mpirun_rsh -np 32 -hostfile $PBS_NODEFILE montecarloparallel

Exercises

1. Explain why using Monte Carlo to evaluate ${\displaystyle \int _{0}^{1}{\frac {1}{1+x^{2}}}\mathrm {d} x}$ allows you to find ${\displaystyle \pi }$ and, in your own words, explain what the serial and parallel programs do.
2. Find the time it takes to run the Parallel Monte Carlo program on 32, 64, 128, 256 and 512 cores.
3. Use a parallel Monte Carlo integration program to evaluate ${\displaystyle \iint x^{2}+y^{6}+\exp(xy)\cos(y\exp(x))\mathrm {d} A}$ over the unit circle.
4. Use a parallel Monte Carlo integration program to approximate the volume of the ellipsoid ${\displaystyle {\frac {x^{2}}{9}}+{\frac {y^{2}}{4}}+{\frac {z^{2}}{1}}=1}$. Use either OpenMP or MPI.
5. Write parallel programs to find the volume of the 4 dimensional sphere ${\displaystyle 1\geq \sum _{i=1}^{4}x_{i}^{2}.}$ Try both Monte Carlo and Riemann sum techniques. Use either OpenMP or MPI.