Parallel Computing in Mathematica

Richard Gass

Department of Physics
University of Cincinnati

Introduction

Parallel Computing in Mathematica is

• 	Very Simple to Configure
•	Once configured it is easy to use
•	Makes it easy to create ad-hoc clusters
•	Much easier than say MPI

However parallel is still harder than serial. Much harder if the problem is not embarrassingly parallel.

Of course not everything can be done in parallel.

Parallel

Most Linear Algebra
FFT's
Monte Carlo Calculations
Numerical Solutions of PDE's

Inherently Serial

Simpification
Any recursive calculation

How much can we gain?

Amdahl's Law

S[p_,f_]=1/(1-f+f/p);

Manipulate[Plot[S[n,f],{f,0,1},FrameLabel->{"fraction of code run in parallel","Speedup"},PlotLabel->"Amdahl's Law",RotateLabel->True,Frame->True,PlotRange->{1,16}],{{n,2,"Processors"},2,16,2,Appearance->"Labeled"},SaveDefinitions->True]

How much do we gain?

Communication Overhead

Things are even worse than Amdahl's law implies due to communication overhead.

Time to send $KernelID and receive the result

Grids and Clusters

• Grids

Loosely coupled collection of machines, may not be homogeneous or in the same geographic location

• Clusters

Tightly coupled, homogeneous, maybe runing managment and queuing software

Set-Up and Kernel Status

Trivial if Lightweight Grid manager is installed.

Set-Up and Kernel Status

Embarrassingly Parallel Examples

Random Number Generation - Compute 16 million random numbers

On a single core

t1=AbsoluteTime[];
data=Table[RandomReal[],{i,1,16 (10^6)}];
timeLocal=AbsoluteTime[]-t1

On four cores

ParallelEvaluate[$KernelID]

t1=AbsoluteTime[];
data=ParallelTable[RandomReal[],{j,1,16},{i,1,(10^6)}];
timeParallel=AbsoluteTime[]-t1

Resulting in a slow-down by a factor of more than 2.6.
Trival computation but lots of communication overhead

Embarrassingly Parallel Examples

Calculation of hard integrals

f[{i_,j_}]:=NIntegrate[(1/25) ((Sin[(1/10) π i (x+5)])^2) ((Sin[(1/10) π j (y+5)])^2) (2.5` Boole[(x^2)+Sin[x y]+((y^2)/4)>1]+(1/2) ((x^2)+(y^2))),{x,-10/2,10/2},{y,-10/2,10/2},PrecisionGoal->3,Method->"AdaptiveMonteCarlo"]

DistributeDefinitions[f]
data=Table[{i,j},{i,1,8},{j,1,8}]
result=WaitAll[Map[Composition[ParallelSubmit,f],Flatten[data,1]]]

Domain Decomposition

Solve ((∂U(t,x))/(∂t))=κ (∇^2)U(t,x)

•	Use Crank-Nicolson to convert the PDE into a set of linear equations

 -r u((Subscript[t, j+1]),(Subscript[x, i+1]))+(2+2r) u((Subscript[t, j+1]),(Subscript[x, i]))-r u( (Subscript[t, j+1]),(Subscript[x, i-1]))  =
r u( (Subscript[t, j]),(Subscript[x, i+1]))+ (2-2r) u((Subscript[t, j]),(Subscript[x, i]))+r u((Subscript[t, j]),(Subscript[x, i-1]))
 r=κ Δt/(((Δx))^2)

•	Yields a tridiagonal system of equations for which LinearSolve is very efficient. 
Keep all vectors and matrices as packed arrays and/or sparse arrays. 
SetSystemOptions["PackedArrayOptions" -> "UnpackMessage" -> True] is useful.

•	Break the spatial domain into pieces. Each processor gets a piece
•	Need communication at the edges
•	Asynchronous method - each processor uses the result from the others one time step back

Hypothetical 30 element grid on three processors

Code

InitialConcentration=Table[((N[Sin[2Pi i/40000]])^2),{i,0,40000}];Clear[u,r,vectors,matrices]
leftend=0.0;
rightend=0.0;
Table[u[0,i]=InitialConcentration[[i]],{i,1,Length[InitialConcentration]}];
u[i_,0]=leftend;
u[i_,Length[InitialConcentration]+1]=rightend;
var[j_]:=Table[u[j,i],{i,Length[InitialConcentration]}]
eqn[j_,i_]:=- u[j+1,i-1]+(2/r+2)u[j+1,i]-u[j+1,i+1]== u[j,i-1]+ u[j,i+1]+(2/r-2) u[j,i]
equations[j_]:=Table[eqn[j,i],{i,1,Length[InitialConcentration]}]
n=Length[ParallelEvaluate[$KernelID]];

width=Ceiling[Length[InitialConcentration]/n];
edges=Table[{i+1,Min[i+width,Length[InitialConcentration]]},{i,0,Length[InitialConcentration],width}];
segments[i_]:=Transpose[{Map[Take[equations[i-1],#]&,edges],Map[Take[var[i],#]&,edges]}]

r=10;arrays=Map[CoefficientArrays[#1[[1]],#[[2]]]&,segments[1]];
vectors=-Map[First,arrays];
matrices=Map[Last,arrays];
DistributeDefinitions[matrices,vectors]
t1=AbsoluteTime[];
Do[u[i,x_]=u[i-1,x];
newU=WaitAll[MapThread[ParallelSubmit[LinearSolve[#1,#2]]&,{matrices,vectors}]];
edgeValues=Map[Drop[#,{2,-2}]&,newU];
boundaries=ReplacePart[edgeValues,{{1,1}->leftend,{n,2}->rightend}];
vectors=MapThread[Take[RotateRight[Insert[Insert[#1,2#2[[1]],1], 2#2[[2]],-1]]+RotateLeft[Insert[Insert[#1, 2#2[[1]],1],2 #2[[2]],-1]],{2,-2}]+(2/r-2) #&,{newU,boundaries}],{i,1,4000}]
t2=AbsoluteTime[]-t1

Errors

Four  processors with asynchronous updating on a 40 point grid. Dashed line is solution on a single processor

Error on a 40,000 point grid.

Performance

Speed-up on Glenn at OSC. Chips are dual socket quad core Opterons