10 YEARS OF MATHEMATICA

Parallel Programming

Roman E. Maeder

MathConsult Dr. R. Mäder
Samstagernstrasse 58a
8832 Wollerau
Switzerland
mailto:maeder@mathconsult.ch
http://www.mathconsult.ch/

The MathLink communication protocol can be used to control several Mathematica kernel processes from within Mathematica. This feature allows the implementation of an environment for parallel programming using processes with distributed memory. Parallel language constructs, such as a parallel version of Map, can easily be implemented on top of these primitive operations.

We describe a library for parallel programming in Mathematica and give several examples of its use, including methods to synchronize global variable updates.

A presentation given at the 1998 Worldwide Mathematica Conference, Chicago, June 1998.

© 1998 MathConsult Dr. R. Mäder. All rights reserved.

Introduction

Parallel computing:

In this context:

Please look for this notebook and the packages
Parallel.m
VirtualShared.m
on my web site, http://www.mathconsult.ch/math/stuff/

MathLink Primitives

Controlling Kernels

Don't evaluate initialization at this point

Start another local kernel (the command to do so may depend on the operating system):

local = LinkLaunch["math -mathlink"]

Start a remote kernel with rsh or ssh:

remote = 
      LinkLaunch[
        "ssh -e none byblis.wolfram.com /home/rd2/scripts/math -mathlink"]

You can also start the remote kernel manually with math -mathlink -linkcreate and then connect to the resulting port using LinkConnect.

Swallowing initial line-number packet.

LinkRead[local]

LinkRead[remote]

Avoid blocking:

LinkReadyQ[local]

Sending expressions for evaluation:

LinkWrite[local, Unevaluated[EvaluatePacket[1 + 2]]]

Polling for result

If[LinkReadyQ[local], LinkRead[local], Null]

Dealing with other packets

LinkWrite[remote, Unevaluated[EvaluatePacket[1/0]]]

If[LinkReadyQ[remote], LinkRead[remote], Null]

Ending session

LinkClose/@{local, remote}

Lining up Workers

Auxiliary routines for starting/closing MathLink connections (OS dependent).

LaunchSlave["host", "cmd_template", options…]

ConnectSlave["port@host", options…]

cmd_template is a string used in StringForm, `1` will be replaced by the hostname. The second form connects to a waiting kernel; this form can be used if security prevents you from executing a remote command directly.

Predefined templates:

$rsh = "rsh `1` math -mathlink"
$ssh = "ssh -e none `1` math -mathlink"

Start local slave (the hostname is ignored):

LaunchSlave["whatever", "math -mathlink"]

List of all slaves: $Slaves

Sending and Receiving Expressions

SetAttributes[Send, HoldRest]

Send[link_LinkObject, expr_]:=
  LinkWrite[ link, Unevaluated[EvaluatePacket[expr]] ]

Options[Receive] = {Wait->False}

Receive[link_LinkObject, opts___?OptionQ] :=
  Module[{wait, res},
    wait = Wait /. {opts} /. Options[Receive];
    While[ wait || LinkReadyQ[link],
      res = LinkRead[link];
      If[ res === $Failed || Head[res] === LinkRead,
            $Slaves = Complement[$Slaves, link];
            Return[$Failed]
      ];
      If[ Head[res] === ReturnPacket,
          Return[res[[1]]] ];
      (* other packets *)
      Handle[link, res];
    ];
    (* nothing there *)
    $NotReady
  ]

Basic remote evaluation (blocking):

SetAttributes[Eval, HoldRest]

Eval[link_LinkObject, expr_]:=
    Module[{res},
        Send[link, expr];
        Receive[link, Wait->True]
    ]

Process Queues

Queue jobs, wait for (some of them) to finish, get results.

Starting Jobs (Processes), Waiting for Results

Queue[cmd]

Enters the given expression into queue, returns a process id (cookie).

WaitOne[pid_list]

Waits for one of the given processes to finish and returns {result,pid,remaining_ids}.

Wait[pid_list], Wait[pid]

Wait for all processes or a single process, return list of results and single result, respectively.

The associative array Id[pid] keeps track of the state of all processes in the various states; first, a process is queued, then it is running on a certain link, then it is finished and waits for collection of result.

queued, link, finished[result]

Enter a job into queue

SetAttributes[Queue, HoldAll]

Queue[cmd_] :=
  Module[{pid},
    AbortProtect[
      AppendTo[$Queue, Hold[cmd, pid]];
      Id[pid] = queued ];
      trace[1, pid, " queued (", Length[$Queue], ")"];
    QueueRun[]
    pid
  ]

Get results of finished jobs:

DoneQ[id_] := Head[Id[id]]===finished

DeQueue[id_] :=
  With[{res=Id[id][[1]]},
    Unset[Id[id]];
    trace[1, id, " dequeued"];
    res
  ]

Wait[id_Symbol] := (
    While[ !DoneQ[id], Pause[$pause]; QueueRun[] ];
    DeQueue[id]
  )

Wait[ids_List] := Wait /@ ids

WaitOne[ids_List/;Length[ids]>0] :=
  Module[{i=1, n=Length[ids]},
    QueueRun[];
    While[True,
      With[{id = ids[[i]]},
        If[DoneQ[id],
           Return[{DeQueue[id], id, Drop[ids,{i}]}],
         (* else check next one *)
           i = Mod[i, n] + 1;
           If[i == 1, Pause[$pause]; QueueRun[]];
        ]
      ]
    ]
  ]

Managing the Queue, Dispatching Jobs

When link becomes available, start a queued process.
When result is waiting, store it and put link back on pool of idle links.

Links are in one of two queues: $Busy, $Idle

Go once through both link queues:

QueueRun[] :=
  Module[{res, i=1},
    (* find finished jobs *)
    While[i <= Length[$Busy],
      Switch[res = Receive[$Busy[[i]], Wait->False],
        {_,_}, AbortProtect[ (* finished *)
               Id[res[[2]]] = finished[res[[1]]];
               trace[1, res[[2]], " done"];
               AppendTo[$Idle, $Busy[[i]]];
               $Busy = Drop[$Busy, {i}] ],
        $NotReady, (* still working: try next *)
               i++,
        _, (* problem, drop link *)
               $Busy = Drop[$Busy, {i}]
      ]
    ];
    (* start queued jobs *)
    While[Length[$Queue]>0 && Length[$Idle]>0,
      AbortProtect[With[{link=First[$Idle]},
        Replace[ First[$Queue],
            Hold[cmd_,id_] :>
              (Id[id]=link; Send[link,{cmd,id}]) ];
        trace[1, First[$Queue][[2]], " started"];
        $Idle=Rest[$Idle];
        AppendTo[$Busy, link];
        $Queue = Rest[$Queue];
      ]]
    ];
  ]

Examples

Needs["Parallel`"]

Create to slaves on local host

slaves = 
      Table[LaunchSlave["localhost", "math -mathlink"], {2}]

Queue a job and wait for it to finish.

j1 = Queue[Unevaluated[1 + 1]]

Wait[j1]

Queue several jobs

Queue/@Hold[1 + 1, 2 + 2, 3 + 3]

pids = List@@%

Current queue states:

{$Idle, $Busy, $Queue}

nondeterministic waiting for a job:

WaitOne[pids]

A convenient use of the return value to drain queue:

pids = %[[3]];

{res, pid, pids} = WaitOne[pids]

Cleanup:

CloseSlaves[]

Quit[];

Parallel Evaluation and Mapping

Evaluating Elements of an Expression in Parallel

ParEval[h[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]]

Sends the (Subscript[e, i]) to available links, then collects results.

We developed the basic idea in the previous example.

SetAttributes[ParEval, HoldAll]

ParEval[h_[exprs___]] :=
      h @@ Wait[ List @@ Queue /@ Hold[exprs] ]

Parallel Mapping

ParMap[f, 
      h[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]]

Evaluates Map[f, h[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]] in parallel.

We could reduce it to ParEval, but it's so simple to express directly:

ParMap[f_, h_[exprs___]] :=
      h @@ Wait[ List @@ Queue /@ f /@ Hold[exprs] ]

Example and Traces

Evaluate initialization here.

links = start[3]

\(ParMap[#\^2&, {1, 2, 3, 4}]\)

$traceLevel = 1;

$traceLevel = 0;

Initializing Slaves

Move definitions of global functions and variables to the slaves:

ExportEnvironment[{(Subscript[f, 1]), (Subscript[f, 2]), …, (Subscript[f, n])}]

ExportEnvironment["context`"]

Mathematica code is portable: no need to have machine-specific code on slaves!
Use existing MathLink connection to transport code: no distribution problems, no common FS needed.

Example: Eigenvalues of Matrices

A parameter and two functions.

prec = 18;

mat[n_] := Table[Random[Real, {-1, 1}, prec], {n}, {n}]

tf[n_] := Timing[Eigenvalues[mat[n]]][[1]]

Make them known to all slaves.

ExportEnvironment[Unevaluated[{prec, tf, mat}]];

Check:

Eval[#, prec]& /@ $Slaves

A sample run

ParMap[tf, {5, 10, 15, 20}]

$traceLevel = 1;

$traceLevel = 0;

Example: an Infinite Calculation

poly[n_, x_] := (∑_(i = 1)^(n + 1))i (x^(i - 1))

ExportEnvironment[{poly}];

d = 50; ids = {};
While[True,
\(While[Length[$Queue] == 0, 
	
                With[{d = d}, 
                  AppendTo[ids, Queue[{d, Length[FactorList[poly[d, x]]]}]]
                    
]; 
                
\(d++\)
]\);
If[Length[ids] > 0, 
                
{res, id, ids} = WaitOne[ids]; 
                
Print[res]
];
	]

Collect orphans after an interrupt:

Wait[ids]

Performance Measurements

Overhead

Time for servicing one trivial process.

With[{n = 200, t = AbsoluteTime[]}, 
      
ParMap[Identity, Range[n]]; 
      
((AbsoluteTime[] - t)/n)]

Speedup and Efficiency

Speedup is the ratio of total parallel execution time to sequential time.

Efficiency is the ratio of observed speedup to number of processors.

Performance[ntry_, d_:  1]:=
Module[{p = Length[$Slaves], t, su},
t = AbsoluteTime[];
ParMap[Pause[d]&,Range[ntry]];
su = ((ntry d)/(AbsoluteTime[] - t));
{p, su, su/p}]

Performance[10]

Add more processors

start[3]; Length[$Slaves]

Performance[24]

Performance[25]

CloseSlaves[]

Quit[];

Virtual Shared Memory

How can processes communicate with each other?

Virtual Shared Memory: faking it on distributed memory machines: PVM.

When a process needs the value of a global (shared) variable or wants to set its value, it sends a request back to the scheduler. The scheduler maintains the variables, udpates their values, and sends back the requested value.

Implementation: send appropriate packets and handle them in the code of Receive, where we already put a handler in place:

Receive[link_LinkObject, opts___?OptionQ] :=
  Module[{wait, res},
    wait = Wait /. {opts} /. Options[Receive];
    While[ wait || LinkReadyQ[link],
      res = LinkRead[link];
      If[ Head[res] === ReturnPacket,
          Return[res[[1]]] ];
      (* other packets *)
          Handle[link, res];
    ];
    (* nothing there *)
    $NotReady
  ]

Slave Side

A shared variable is a symbol s with two definitions attached to it, one for reading the value, one for setting a new value (an upvalue on Set).

s := `Slave`GetValue[s];
s/: (s=e_) := `Slave`GetValue[s, e];

Similarly, shared downvalues (associative arrays, functions), have appropriate definitions:

d_s := `Slave`GetValue[d];
s/: (d_s=e_) := `Slave`GetValue[d, e];

Note that only ordinary definitions with = are trapped. Delayed values and conditional rules  introduce scoping nightmares.

The auxiliary function GetValue performs the callback and waits for the requested value. The request is wrapped in CallBackPacket, the reply arrives wrapped in ReplyPacket.

SetAttributes[{`Slave`GetValue, `Slave`CallBackPacket}, HoldFirst];

c_`Slave`GetValue :=
  Module[{r},
    LinkWrite[$ParentLink, `Slave`CallBackPacket[c]];
    While[Head[r = LinkRead[$ParentLink]] =!=
                   `Slave`ReplyPacket, Null];
    r[[1]]
  ]

Master Side

The hook Handle is overloaded to deal with CallBackPacket. The packet contains the original GetValue request; one argument is a lookup, two arguments a variable set.

`Slave`CallBackPacket/:
Handle[link_, `Slave`CallBackPacket[`Slave`GetValue[s_,e_:ReadOnly,t_:False]]] :=
  Module[{},
    Which[
      t, (* test and set *)
          If[!ValueQ[s] || s === Null,
               s = e;
               trace[2, HoldForm[s ? e]],
             (* else: locked *)
               trace[2, HoldForm[s ? e], " (", s, ")"]
              ],
      e =!= ReadOnly,  (* set value *)
          s = e;
          trace[2, HoldForm[s = e]],
      True, (* get value *)
          trace[2, HoldForm[s]->s]
    ];
    If[!ValueQ[s], s = NoValue];
    LinkWrite[link, `Slave`ReplyPacket[s]];
  ]

To avoid infinite recursion, all variables need a value. Unitialized variables are set to NoValue.

Declaring Shared Variables/Functions

A declaration sets up the slave-side definitions and transports them to all slaves. The definition for GetValue is also assembled and shipped to the slaves, which, therefore, need no separate init code.

SharedVariables[{(Subscript[v, 1]), (Subscript[v, 2]), …, (Subscript[v, n])}]

SharedDownValues[{(Subscript[f, 1]), (Subscript[f, 2]), …, (Subscript[f, n])}]

Examples

Re-init kernel here...

<< VirtualShared`

links = start[3]

Declare shared variables.

SharedVariables[{x, y, z}]

You can see the slave-side definitions:

Eval[links[[1]], Information[x]]

Basic tests.

x = 0;

Eval[links[[1]], x]

Eval[links[[1]], x = x + 1]

x

Eval[links[[2]], x]

The prototypical synchronization problem: between reading and writing, another process can change the value or read the obsolete value.

y = 0;
ParMap[(
		Pause[0.5 Random[]];
		a=y;
		Pause[Random[]];
		y=a + 1)&,
	Range[10]]

y

This problem appears also on multiprocessing or timeshare operating systems.

Synchronization

Exclusive access to variables can be achieved with one new primitive operation, a conditional atomic write. In machine code, it is often called test and set or read-modify-write.

TestAndSet[s, e]

Set s to e (just like s=e), but only if s is currently unset or its value is Null. Return the (old or new) value of s. The read-modify-write is atomic.

The slave-side definition is similar to the one for Set (=):

s/: Global`TestAndSet[s, e_] :=
        `Slave`GetValue[s, e, True];
s/: Global`TestAndSet[d_s, e_] :=
        `Slave`GetValue[d, e, True];

The additional code in the handler was given earlier.

links = $Slaves

SharedVariables[{z}]

Eval[links[[1]], TestAndSet[z, 17]]

Eval[links[[2]], TestAndSet[z, 18]]

Eval[links[[1]], z = Null]

Eval[links[[2]], TestAndSet[z, 18]]

Critical Sections

Exclusive access to a variable or other resource can be controlled by a lock. Here is the typical code added to the previous example about asynchronous access to a variable.

SharedVariables[{y, lock}]

$traceLevel = 2;

y = 0;
ParMap[(
		Pause[0.5  Random[]];
		While[TestAndSet[lock, #] != #, Pause[0.2]];
		a = y;
		Pause[Random[]];
		y = a + 1;
		lock = Null;
		a + 1)&,
	Range[10]]

CloseSlaves[]

The Dining Philosophers

The classic example for access to shared resources, in this case forks placed between philosophers sitting at a round table. A philosopher needs the two forks to his left and right to eat.

start[4]

startRemote[4]

The variables forks[1],…,forks[n] serve to reserve forks by conditionally setting them to a unique value.

SharedDownValues[{forks}]

stop = False; SharedVariables[{stop}]

Each dining philosopher is one process.

n = Length[$Slaves]

Here's the behaviour of a philosopher: either thinking or eating.

philo[i_]:=
With[{f1 = i, f2 = Mod[i, n] + 1},
While[! stop,
Print[i, " thinking"];
  
Pause[10 Random[]];
		  
Print[i, " hungry"];
			While[TestAndSet[forks[f1], i] != i || 
                        TestAndSet[forks[f2], i] != i,
Print[i, " angry"];
Pause[0.3]];
  
Print[i, " eating"];
Pause[3 Random[]];
  
forks[f1] = Null;forks[f2] = Null;
]]

ExportEnvironment[Unevaluated[{philo, n}]];

We want essentially ParMap, but need a way to stop them cleanly.

ParMap[philo, Range[n]]

pids = List@@Queue/@philo/@Hold@@Range[n]; Wait[pids]

stop = True; Wait[pids]

CloseSlaves[]

Latency Hiding

Latency: communication overhead, the wait between end of servicing a request and start of servicing new one.

MathLink is a buffered stream. Can send additional requests before receiving all outstanding ones.
Keep buffers full at all times: no latency.

Does not work with current implementation of virtual shared variables.

Latency Hiding is not implemented in the demo version of the Parallel.m packages available on my Web server.

Example, Revisited

Go back to polynomial factoring with

$loadFactor = 2;

start[1]

startRemote[Length[wriHosts]]

Eval[#, {$MachineName, $ProcessID, $Version}]&/@$Slaves // TableForm

Conclusions

Further topics in parallel processing:

Scratch

Start

Init Code

<< Parallel`

Start up Slaves

Local ones

start[n_] := 
      Table[LaunchSlave["localhost", "math -mathlink -noinit"], {n}]

put available hosts here

wriHosts=ToString/@{};

startRemote[n_] /;  n <=  Length[wriHosts] := 
      LaunchSlave[#, $wsh]& /@ Take[wriHosts, n]

n local slaves

links = start[4]

n remote slaves

links = startRemote[4]

All available ones

links = startRemote[Length[wriHosts]]

Status

$Slaves

TableForm[{$Idle, $Busy, $Queue}, TableDepth -> 1]

?? Id

Eval[#, {$MachineName, $ProcessID, $Version}]&/@$Slaves // TableForm

Stop

CloseSlaves[]

Quit[];