Enhancements in Parallel Computing in Mathematica

Roman E. Maeder

Wolfram Research, Inc.

From a 2 core laptop to a supercomputer with 100 thousands of processors, there is a wide variety of hardware available for parallel computing. The challenge is to find a uniform approach that adapts to most configurations.
Such tools for parallel computing are now built into every copy of Mathematica.

The History of Parallel Computing in Mathematica

The Prehistory

Parallel Computer Algebra

MathLink

Parallel Computing Toolkit 1.0

Parallel Computing Toolkit (PCT) was first available in 1997. At that time we had what we would now call 2 core workstations from Sun and otherwise relied on a 10MB/s ethernet for connection to other machines.

PCT's method to acquire remote kernels was to specify an OS command to remote login and start a kernel.

$RemoteCommand

Starting local kernels is straightforward on any OS.

$mathkernel

Features

Symbolic, treats programs as data

Fault-tolerant

Very flexible

Virtual Shared Memory

Parallel Computing Toolkit 2.0 and gridMathematica

Automatic Parallelization

ParallelEvaluate[Map[f,{a,b,c}]]

Deployment

gridMathematica for bulk licensing of batch-operated kernels.

Ease of Use

Cluster Integration Package to interface with cluster management systems.

The Next Version of Mathematica

Integration

Out of the box support for multi-core workstations

Parallel Tools (the former Parallel Computing Toolkit) is integrated into the core language

Context[ParallelMap]

Single-Machine License supports multi-core machines

Ease of Use

ParallelMap[PrimeQ[2^# - 1] &, Range[30000, 30031]]

Parallelize[Map[PrimeQ[2^# - 1] &, Range[30000, 30031]]]

Configuration and Monitoring

GUI-based configuration and monitoring tools.

Additional Kernel Connection Methods supplied

For accessing various computer resources.

Workbench Integration

For developing and debugging parallel programs.

Environments for Parallel Computing

>From a 2 core laptop to a supercomputer with 100 thousands of processors, there is a wide variety of hardware available for parallel computing. The challenge is to find a uniform interface that adapts to most configurations.

Hardware Environments

Multi-Core Computers

This is now standard, and we support it inside standard Mathematica.

$ProcessorCount

Clusters

Managed or unmanaged racks of identical components.

Ad-hoc Networks

Computer labs, networked workstations, remote machines.

Kernel Connection Methods

Local Kernels on multi-core computers

Uses a shared-memory implementation of MathLink for efficiency.
By default launches one parallel kernel per core.

LaunchKernels[4]

Remote Services

Finding available Mathematica resources on ad-hoc networks; resources advertise themselves.

The client is part of Mathematica and interfaces nicely with the configuration tools.

Cluster Integration

Big, managed clusters do not allow just anyone to log into their nodes.

Our solution: interface with the cluster management software to launch kernels and then pass them to Parallel Tools.

The client is part of Mathematica and interfaces nicely with the configuration tools.

Remote Kernels (remote login)

Lets you configure commands and hostnames for automated launching of remote kernels; the standard method in PCT.

Types of Parallel Computing

Doing the same thing many times

A large percentage of what people do with Mathematica involves computations with many independent subtasks.

Performing the same computation with many different inputs

doing a Monte-Carlo simulation

rendering frames of an animation

In a functional style, all of these involve Map[]. To parallelize, simply use ParallelMap[] instead.

Data parallelism

Many operations of the form f[some long list] can be parallelized using this method:

divide the list into as many pieces {(Subscript[e, i]), …, (Subscript[e, j])} as there are processors

compute f[{(Subscript[e, i]), …, (Subscript[e, j])}] in parallel (giving (Subscript[r, k]))

combine the results comb[(Subscript[r, 1]), …, (Subscript[r, m])]

ParallelCombine[f,Range[20],comb]

Concurrency

More general parallelism for less structured problems that requires communication among processes.

Note that a process in this context is simply a Mathematica expression to be evaluated (and a processor is a running Mathematica kernel).

A scheduler sends waiting processes to available processor and collects the results.

Domain Decomposition

For the hard cases: as single problem such as a partial differential equation that requires tight interaction among the processors that work on part of the domain. The realm of MPI.

Parallel Tools

All language constructs to support parallel programming are part of Mathematica.

No other code or common file system is necessary.

The notion of programs as data makes even failure recovery rather easy to implement.

Parallel Tools interfaces nicely with the Workbench, allowing debugging of remote kernels and parallel computations.

Compatibility with Parallel Computing Toolkit

If you ever developed programs for Parallel Computing Toolkit 2.1, remember this command:

Needs["Parallel`V2`"]

It gives you back all the old command names.

Also note that Version 6 subkernels are supported from Version 7 master kernels.

Data Parallelism

A large number of computations involve independent operation on the elements of a long list. This includes the independent subtasks mentioned above, where the list is implicit (such as a list of random numbers).

ParallelTable[$KernelID,{12}]

We handle this class with ParallelCombine to split the long list into pieces.

Divide et Impera

Testing numbers for primality

Parallelize[PrimeQ[{1,2,3,4,5,6,7,8,9,10,11}]]

Parallel pattern matching

Parallelize[Cases[Range[100],_?PrimeQ]]

Membership testing and counting

Parallelize[Count[Range[1000],_?PrimeQ]]

Automatic Parallelization

Parallelize recognizes expressions that are data parallel.

Parallelize[Outer[f,{a,b,c},{x,y,z}]]

Ultimate goal

Parallelize[your program goes here...]

Concurrency

The Life of a Process

pid=ParallelSubmit[1+2]

WaitAll[pid]

Watch the scheduler in action.

pids=ParallelSubmit[PrimeQ[2^#-1]]&/@{1,2,3,4,5,6,7}

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

WaitAll[pids]

SetOptions[$Parallel,Tracers->{}]

Functional Programming and Process Generation

WaitAll[ParallelSubmit[expr]] is just the evaluation function (or identity).

WaitAll[ParallelSubmit[1+2]]

To evaluate a function remotely, you can use

WaitAll[Composition[ParallelSubmit,f][x]]

This observation gives us a powerful recipe for parallelizing functional operations. This example develops a parallel version of the Mathematica function Map[].

Map[f,{a,b,c}]

We compose f with ParallelSubmit[]. The result is a list of evaluations that represent the individual function applications.

Map[Composition[ParallelSubmit,f],{a,b,c,d,e,g}]

WaitAll[] runs the scheduler and returns with the results, which were computed in parallel on any available parallel kernels.

WaitAll[%]

The powerful functional language allows us to derive parallel code automatically.

Speculative Evaluation and Nondeterminism

After submitting a number of processes the master kernel usually has to wait for at least one of them to finish.

It does not need to wait for all of them to finish, however. Processes can be aborted (a kind of parallel ⌘-.)

Searching

Finding an irreducible polynomial by testing random polynomials

Finding a prime in a given range

In general: find an element of an implicit or explicit search space with a desired property

Distribute the search space and run concurrent searches.
As soon as a result comes back, abort all remaining processes.

Sequential Code

Generate a random polynomial of degree n mod p.

randPoly[n_Integer?Positive,x_,p_]:=Random[Integer,{1,p-1}]*x^n+Sum[Random[Integer,{0,p-1}]*x^i,{i,1,n-1}]+Random[Integer,{1,p-1}]

Test for irreducibility.

irreducibleQ[poly_,p_]:=IrreduciblePolynomialQ[poly,Modulus->p]

An example mod 2.

poly=randPoly[10,x,2]

irreducibleQ[poly,2]

Generate new polynomials until one is irreducible.

findOne[n_,p_,x_:x]:=Module[{poly},
While[!irreducibleQ[poly=randPoly[n,x,p],p],Null];
poly]

findOne[100,2]

Distribute Code

Distribute the code to all remote kernels (Mathematica code is data, so it can be collected and sent across MathLink).

DownValues[randPoly]

DistributeDefinitions[randPoly,irreducibleQ,findOne]

Parallel Execution of the Search Loops

Run the search loop on each remote kernel.
As soon as one result arrives we can abort the remaining remote kernels.

ParallelTry[findOne[1000,2]&,Range[$KernelCount]]

Clear[poly]

Try out different methods

Compare different ways of solving the same problem.

With[{a=Reverse/@Distribute[{{-32,-16,0,16,32},{-32,-16,0,16,32}},List]},f[x_,y_]:=1/(0.002+Plus@@MapIndexed[1/(#2[[1]]+Plus@@(({x,y}-#1)^6))&,a])
];
ranges={{x,-50,50},{y,-50,50}};

Plot3D[f@@First/@ranges,Evaluate[Sequence@@ranges]]

This function takes a method as arguments and solves a hard-coded problem using it.

tryMethod[method_]:=
Module[{time,res},
{time,res}=Timing[NMinimize[f[x,y],ranges,Method->method,PrecisionGoal->12]];
Join[{method/.{meth_,___}:>meth,time},res]
]

methods={
"DifferentialEvolution",
"NelderMead",
"RandomSearch",
"SimulatedAnnealing",
Automatic
};

tryMethod[methods[[1]]]

Compare Methods In Parallel

DistributeDefinitions[f,ranges,tryMethod]

TableForm[ParallelMap[tryMethod,methods],TableDepth->2,TableHeadings->{None, {Method,Timing,Minimum,Position}}]

Let the first one win

ParallelTry[tryMethod,Complement[methods,{Automatic}]]

Shared Resources

For non-trivial problems, there is some communication among the running processes.

Some resources are shared (sequential) and access needs to be synchronized.

Distributed Values

Distribute the current definitions of variables and functions to all kernels, but let each one maintain its own independent local copy thereafter.

Used for distributing code and parameters.

DistributeDefinitions[variable, function,…]

ParallelNeeds["Package`"]

Kernels can be initialized solely through the existing MathLink connection, no common file system or code distribution is necessary.

Shared Variables

Variables global to the computation. There is only one copy, maintained by the master kernel. All access is through callbacks to the master kernel.

SetSharedVariable[xs]; xs=0;

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

ParallelEvaluate[xs=xs+1,Kernels[][[1]]]

ParallelEvaluate[xs]

SetOptions[$Parallel,Tracers->None];
UnsetShared[xs];

Synchronization

What if two kernels do x=x+1 at the same time?

SetSharedVariable[xs]

xs=0;
ParallelEvaluate[xs=xs+1]
xs

UnsetShared[xs];

What if two kernels write to the same file at the same time?

Locking

Parallel Computing Toolkit contains all primitives needed to implement locking and exclusive access to critical sections.

Here, we run five processes that all access a shared variable ys, first reading it and some time later updating it. The code in between is protected by obtaining a lock and releasing it afterward.

SetSharedVariable[ys]
Clear[lock]

ys=0;
ParallelMap[Module[{a},
	Pause[Random[]];
	CriticalSection[{lock},
		a=ys;
		Pause[0.5Random[]];
		ys=a+1;
	];
	a+1]&,
  Range[10]]

UnsetShared[ys];

Shared Definitions

Shared variables for downvalues (definitions of functions).

Clear[f]; SetSharedFunction[f]

ParallelEvaluate[f[$ProcessID]=$KernelID]

?f

UnsetShared[f];

Can even do parallel dynamic programming.

Callbacks, Threads

Shared variables are essentially implemented via shared downvalues. The callback is atomic (thread-safe), and can therefore be used for synchronization.

Example: Parallel Sow and Reap

Parallelizing Sow/Reap: sow on the parallel kernels, reap all results in the master kernel. But the calls of Sow[] need to be executed on the master, not the kernels where they are generated. We define an auxiliary shared downvalue sow[e].

sow[e_]:=(Sow[e];Null)
Protect[sow];
SharedDownValues[sow]

The extra Null prevents the value being needlessly sent back to the parallel kernel where it came from.

Reap[
ParallelEvaluate[sow[$ProcessorID]]
]

UnsetShared[sow]

Great for parallel searching.

Example: Access to a Shared Resource

Run an infinite loop and report back results as they are found. The results are added to a single master list of all results.

Callbacks are implemented with shared functions. They are an alternative tool for synchronization.

Example: tally the number of factors of (2^i)-1.

Clear[count];
count[nf_]:=0

report[deg_,nf_]:=(count[nf]=count[nf]+1;)
Protect[report]; SetSharedFunction[report];

maxcount=20;
Dynamic[BarChart[Table[count[i],{i,maxcount}]]]

WaitAll[Table[ParallelSubmit[{deg},report[deg,Total[FactorInteger[2^deg-1][[All,2]]]]],{deg,1,210}]];

Clear[count];
UnsetShared[report];

Performance

Great new toy: Dynamic variants of parallel inspection features.

Live Evaluation Objects

Internal Hacks

Simplify`PiecewiseThread;
PiecewiseDump`$ShowIntegrals=True;

IntegrateWithAssumptions[integrand_,dx__,opts___?OptionQ]:=Module[{ineqs=#2<#1<#3&@@@{dx},conds,simp=Simplification/.{opts}/.Options[IntegrateWithAssumptions],assum=Assumptions/.{opts}/.Options[IntegrateWithAssumptions]},conds=And@@@NestList[Most,Most[ineqs],Length[ineqs]-1];
Fold[simp[Integrate[#1,#2[[1]],Assumptions->assum&&#2[[2]]],assum&&#2[[2]]]&,integrand,Transpose[{Reverse[{dx}],conds}]]]
Options[IntegrateWithAssumptions]={Simplification->Simplify,Assumptions->True};

Setup

The area of a random triangle in the unit square. The integral involves many cases.

Length[pieces=List@@Integrate[Abs[(-x2 y1+x3 y1+x1 y2-x3 y2-x1 y3+x2 y3)],{x1,0,1},{x2,0,1},{x3,0,1},{y1,0,1},{y2,0,1},{y3,0,1}]]

First[pieces]

First[%/.Integrate->IntegrateWithAssumptions]

Run it

$KernelCount

DistributeDefinitions[IntegrateWithAssumptions]

ids=ParallelSubmit@@@(pieces/.Integrate->IntegrateWithAssumptions)

res=WaitAll[ids]

Simplify[Plus@@res]

Speedup and Methods

The Method Option

All higher-level parallel operations, such as ParalleMap, take a Method option that allows to fine-tune the division into concurrent evaluations to the type of problem.

Many Similar Evaluations

Show the effect of sequential overhead.

ParallelEvaluate[ClearSystemCache["Numeric"]];
ParallelMap[(N[Sin[#],1000];)&,Range[1000],Method->"CoarsestGrained"];

ParallelEvaluate[ClearSystemCache["Numeric"]];
ParallelMap[(N[Sin[#],1000];)&,Range[1000],Method->"FinestGrained"];

Uneven Evaluations

Show the effect of load balancing.

nf=Function[i,IrreduciblePolynomialQ[Cyclotomic[50i+1,x]]];
range=Range[16,1,-1];

ParallelMap[nf,range,Method->"FinestGrained"]

ParallelMap[nf,range,Method->"CoarsestGrained"]