Mainstream Parallel Computing

Experiences with Mathematica 7

Roman E. Maeder

Wolfram Research, Inc.

Mathematica 7 included parallel computing as a standard part of the system, thus offering easy access to parallelism for everyone with a multi-core computer.
We look at experiences gained from user feedback and discuss enhancements forthcoming with the next version of Mathematica, as well as ideas for new features.

Parallel Computing in Mathematica

Numerics: MKL

A large (10000x10000) eigenvalue calculation on 12 cores (2 sockets with six-core Opteron 2431 processors).

Parallel Tools

LaunchKernels[8]

KernelStatus[]

Data Parallelism

Solve and plot a differential equation for many initial conditions and animate the results:

eqn=D[z[t,s],t,t]==-((z[t,s])/(((((z[t,s])^2)+((((1/2) (1+0.1 Sin[2 Pi t])))^2)))^(3/2)));
DistributeDefinitions[eqn]

sols=ParallelTable[localsol=NDSolve[{eqn,z[0,s]==param s,Derivative[1,0][z][0,s]==param},z,{t,0,10}, {s,1,2}]; 
Plot3D[Evaluate[z[tp, z0]/.localsol],{tp,0,10},{z0,1,2},PlotRange->{-5,5},PlotLabel->param], 
{param, -1,1,0.1}];

ListAnimate[sols]

Apply different algorithms to the same set of data:

img=;

Distribute the image to all subkernels:

DistributeDefinitions[img];

Apply a list of different filters to the same image and display the result:

filters={Identity,Sharpen, Blur,GaussianFilter[#,3]&,LaplacianFilter[#,1]&,LaplacianGaussianFilter[#,2]&,
GradientFilter[#,2]&,MedianFilter[#,5]&,
ImageConvolve[#,BoxMatrix[1]/9]&};

ParallelTable[Tooltip[filter[img],filter],{filter,filters}]

Tabulate histograms of word lengths in various languages:

languages=DictionaryLookup[All];Length[languages]

hist=ParallelTable[Histogram[StringLength/@DictionaryLookup[{l,All}],{1,25,1},PlotLabel->l,PlotRange->{{0,25},All},PerformanceGoal->"Speed"],{l,languages}];

Show the results in a grid:

GraphicsGrid[Partition[hist,3],ImageSize->500]

Concurrency

Run a search for a random prime on each parallel kernel:

With[{base=10^1000,r=10^10},WaitAll[Table[ParallelSubmit[
While[!PrimeQ[p=RandomInteger[{base,base+r}]],Null];p],{$KernelCount}]]-base]

Watch the scheduling of the evaluations while they are running:

With[{base=10^1000,r=10^10},evs=Table[ParallelSubmit[
While[!PrimeQ[p=RandomInteger[{base,base+r}]],Null];p-base],{2*$KernelCount}];
PrintTemporary[evs];
WaitAll[evs]]

Nondeterminism

Search for a Mersenne prime starting at a given prime exponent:

search[n0_,dn_]:=
Module[{n=n0},While[!PrimeQ[2^Prime[n]-1],n+=dn];n]
DistributeDefinitions[search]

With[{dn=$KernelCount,n0=500},
n=ParallelTry[search[n0+#,dn]&,Range[0,dn-1]];
{n,Prime[n]}
]

Misunderstandings and Pitfalls

There are many concepts that do not appear in sequential computation.

Communication: Trivial Examples

AbsoluteTiming[par=ParallelTable[Sin[x],{x,0.0,3.14,0.01}];]

AbsoluteTiming[seq=Table[Sin[x],{x,0.0,3.14,0.01}];]

Non-Trivial Example

ParallelTable[StreamPlot[ {(x^i) (y^j),(x^j)(y^i)},{x,-3,3},{y,-3,3},ImageSize->100],{i,1,3},{j,0,3}]

Uneven Loads

If all expressions scheduled for concurrent evaluation will take the same amount of time, creating large batches to send to kernels is best. However, if the times are not all the same, submitting many more smaller batches is better, because it achieves load balancing. Kernels that finish a work load early get assigned a new work load.

Parallelize[Map[Composition[Framed,FactorInteger],{1,11,111,1111,11111,111111}]]

Load Balancing vs. Communications Overhead

Choose the communications pattern based on the property of your calculation.

Break the computation into the smallest possible subunits:

Parallelize[Map[Labeled[Framed[#],$KernelID]&,Range[16]],Method->"FinestGrained"]

Break the computation into as many pieces as there are available kernels:

Parallelize[Map[Labeled[Framed[#],$KernelID]&,Range[16]],Method->"CoarsestGrained"]

Middle ground:

Parallelize[Map[Labeled[Framed[#],$KernelID]&,Range[16]]]

Two Examples

Calculate and display the Feigenbaum (or bifurcation) diagram of the logistics map:

line[r_,dy_,np_,n0_,n_]:=Module[{pts},
With[{logistics=Function[x,r x(1-x)]},
pts=Join@@NestList[logistics,Nest[logistics,RandomReal[{0,1},np],n0],n-1]];
Log[1.0+BinCounts[pts,{0,1,dy}]]]
DistributeDefinitions[line]

With[{w=600,h=350,r0=2.95,r1=4.0},
ArrayPlot[ParallelTable[line[r,1/(w-1),w,1000,100],{r,r0,r1,(r1-r0)/(h-1)}],ImageSize->{w,h},PixelConstrained->True]]

Filled Julia plots (contributed by Roger Germundsson)

ClearAll[ParallelFilledJuliaPlot];
Options[ParallelFilledJuliaPlot]=Sort[Join[Options[Graphics],{PlotPoints:>500,MaxIterations->12}]];
ParallelFilledJuliaPlot[e_,{z_,zmin_,zmax_},o:OptionsPattern[]]:=
Module[{pp,nu,nv,s,u,v},
(* plot points *)
pp=OptionValue[PlotPoints]-1;
If[IntegerQ[pp],{nu,nv}={pp,pp},{nu,nv}=pp];
(* range *)
With[{umin=Re[zmin],umax=Re[zmax],vmin=Im[zmin],vmax=Im[zmax]},
(* iterations *)
With[{du=(umax-umin)/nu,dv=(vmax-vmin)/nv,mi=OptionValue[MaxIterations]},
s=ParallelTable[Log[0.1+Abs[NestWhile[Function[z,e],u+I v,Abs[#]<2&,1,mi]]],{v,vmin,vmax,dv},{u,umin,umax,du}]];
(* display *)
ReliefPlot[s,DataRange->{{umin,umax},{vmin,vmax}},FilterRules[{o},Options[Graphics]],FrameTicks->Automatic,PlotLabel->e]
]]

ParallelFilledJuliaPlot[(z^3)+0.59+I 0.43,{z,-2-2I,2+2I},PlotPoints->800,ImageSize->400]

Manipulate[ParallelFilledJuliaPlot[z^n+Complex@@l,{z,-r(1+I),r(1+I)},ImageSize->300,PlotPoints:>ControlActive[75,300],MaxIterations->i,PlotLabel->{z^n+Complex@@l,i}],
{{n,2},{1,2,3,4,5,6},ControlType->SetterBar},
{{i,10},1,50,1},
{{r,2},1,2},
{{l,{0.,1.}},{-2,-2},{2,2},Locator}]

Undefined Functions

FactorInteger[111111111]

factors[n_]:=Total[FactorInteger[(10^n-1)/9][[All,2]]]

factors[9]

AbsoluteTiming[ParallelTable[factors[n],{n,35,65}]]

DistributeDefinitions[factors]

AbsoluteTiming[ParallelTable[factors[n],{n,35,65},Method->"FinestGrained"]]

Side Effects

list=Table[0,{20}];DistributeDefinitions[list]

ParallelDo[list[[i]]=i^2,{i,Length[list]}]

list

may I suggest:

list=ParallelTable[i^2,{i,20}]

Shared Variables

If you need side effects and communication.

list=Table[0,{20}];
SetSharedVariable[list]

ParallelDo[list[[i]]=i^2,{i,Length[list]}]

list

but this example slows down to essentially sequentical evaluation.

Result Communication

Watch the results appear as they are found:

SetSharedVariable[primes]

Monitor[primes={};
ParallelMap[If[PrimeQ[2^#-1],AppendTo[primes,#]]&,Range[1000,6000]];primes,
primes]

For an alternative way to access a shared resource (such as the variable primes above), see the next section.

Performance & Communication

Real-time Performance Data

Use a shared function to record timing results as they are generated:

record[k_,t_]:=AppendTo[counts[[k]],t];
SetSharedFunction[record]

Set up a dynamic bar chart with the timing results:

ids=ParallelEvaluate[$KernelID];mi=Max[ids];
counts=Table[{0},{mi}];
Dynamic[BarChart[counts[[ids]],PerformanceGoal->"Speed",ChartLayout->"Stacked",ChartLabels->{ids,None}]]

Run a series of calculations with vastly varying runtimes:

counts=Table[{0},{mi}];
Parallelize[Table[Module[{t},t=Timing[Plus@@FactorInteger[2^i-1][[All,2]]];
record[$KernelID,t[[1]]];t[[2]]],
{i,212,180,-1}],Method->"FinestGrained"]

Input and Output Size

Communication for input and output can vary greatly according to subtle variations in problems.

Long Input and Output

ParallelMap[fun,{a,b,c,d,e,f,g,h,i,j}]

Short Output

Parallelize[1+2+3+4+5+6+7+8+9]

Short Input

ParallelTable[fun[i],{i,32},Method->"CoarsestGrained"]

SetOptions[$Parallel,Tracers->{SendReceive}];

ParallelTable[fun[i],{i,32},Method->"CoarsestGrained"];

SetOptions[$Parallel,Tracers->{}];

The trick here is to send a program to the subkernels, rather than expand the table beforehand. Not many systems can do that.

Short Input and Output

SetOptions[$Parallel,Tracers->{SendReceive}];

ParallelSum[i,{i,0.0,1.0,0.00001},Method->"CoarsestGrained"]

SetOptions[$Parallel,Tracers->{}];

Code Analysis (coming soon...)

Hidden tables and ranges:

ParallelMap[Sin,Range[0.0,1.0,0.02]]

Let’s see what it does with it:

SetOptions[$Parallel,Tracers->{SendReceive}];

ParallelMap[Sin,Range[0.0,1.0,0.02],Method->"CoarsestGrained"];

SetOptions[$Parallel,Tracers->{}];

Large Calculations

Unattended Operation

Invoke the master kernel (without frontend) from a script, launch subkernels, and read your program.

You can develop your program in a Mathematica notebook, then automatically save it as a program text file (.m), and load that into a standalone kernel.

Fault Tolerance

Parallel tools supports failure recovery. If a kernel fails for any reason (internal or external), the evaluation running on it is reassigned to another kernel, and the kernel is restarted.

A graph of the memory use on a rather large, 16GB, 12 core, workstation, running a continous parallel Mathematica program. Increased memory use over time is unavoidable here, so eventually, the 12 kernels will hit a system limit and quit.

But they are immediately restarted!

So far, the total elapsed time is 7 days, representing 84 days of CPU time.

Setup

Needs["Parallel`Developer`"]

$ProcessorCount

LaunchKernels[8]

UnsetShared["*"]

Quit[]