Introduction To Parallel Computing

Richard Gass
Department of Physics, University of Cincinnati
Cincinnati, OH 45221-0011
gass@physics.uc.edu

Note

The code in this notebook with work with version 7.01 or 8.0 of Mathematica. Most the code will run in earlier versions but configuring the kernels will be harder.

Why compute in parallel?

In the old days (before 2004) clock speed doubled every 18 months or so and computer manufacturers pushed what is known as the megahertz myth, the faster the clock speed the faster the computer. Of course this is not true. All sorts of things effect performance. The depth of the instruction pipeline, cache size, number of operations per clock cycle and the number of cores all effect performance. Now a days clock speeds are no longer increasing but most computers have more than one core per processor and many have more than one processor. This is being written on a machine with two dual core 2.66 Intel Xenon processors so I have a total of four cores. The Ohio Supercomputer Center's new supercomputer Glenn has over 10,000 cores. Each core on Glenn is about the same speed as the one's on my machine but I only have four. 

If I just want to browse the web, use word, send email and maybe run a single Mathematica kernel then I can let my operating system thread processes over my processors so that the load is balanced and I don't have to think about it. But if I want to use Glenn, a cluster or even just run several Mathematica kernels I have to do parallel computing.

What can we compute in parallel?

Unfortunately not everything can be computed in parallel. A classical example is simplifying a large algebraic expression. There is no known parallel algorithm for this. Other examples include recursive calculations since each step depends on previous steps there is no way to do the calculation in parallel. Thus we can not speed up the calculation of the Fibonacci numbers by doing the calculation in parallel.

fibonacci[n_Integer]:=fibonacci[n]=fibonacci[n-1]+fibonacci[n-2]
fibonacci[2]=1;
fibonacci[1]=1;

fibonacci[200]

However many common operations can be done in parallel. Most linear algebra operations can be done in parallel as can FFT's, most image processing manipulations, Monte Carlo calculations and the numerical solution of many partial differential equations. 

Another type of parallelism called data parallelism occurs when you have many large chunks of data all of which you want to do the same thing to.  In this case you can just ship each chunk off to a different processor.  Explorations in parameter space also lend them selves to parallel computation.

How much do we gain?

If you run your code on n processors you will never see a speed-up by a factor of n although in some cases you may come close. Several factors determine the speed-up the most important being the fraction of code that can be run in parallel and the communication overhead. If f is the faction of code that can be run in parallel and p the number of processors then the potential speed-up is given by Amdahl's law (see Landau, Páez and Bordeianu [1] for a derivation) and is

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

By factor of code I mean of course the fraction of time your code spends doing operations that can be done in parallel.  As you can see from the Manipulate below the results are not encouraging. If you can write 75% of your code in parallel and you run on 16 processors you only get a speed up of just over a factor of three and going to 64 processors buys you very little over 16 processors.

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]

So unless you code spends most of its time doing operations that can be run in parallel you are out of luck. The curve is quite steep for if f=.9 and you run on 16 processors you can expect a speed-up of move than a factor of six.

In practice the situation is worse than Amdahl's law predicts because we have neglected the communications overhead. The processors have to communicate with each other and this takes time. How much time depends both on how much information needs to be exchanged and on what type of cluster or grid or machine you are running on. In the best case, a set of calculations or data are shipped off to the processors by the controlling processor and no communication is needed except to send the results back to the control processor. In cases with high communication overhead compared to computation time per process more processors can actually lead to a slowdown rather than a speedup.

Different types of "clusters and grids"

The terms cluster and grid are often used interchangeably but technically clusters are tightly coupled machines or processors and grids are loosely coupled, not necessarily homogeneous and may be geographically dispersed. Grids may also be collections of clusters as for example TeraGrid is.

Single machines

A multiprocessor single machine is the simplest example of a cluster. The machine I an writing this on has two Duel-Core Intel Xenon chips so I have a total of four cores. That's not very many cores but they are convenient and an excellent tool for debugging parallel code. They also have very low communication overhead since all the cores sit on the same motherboard.

I current have eight kernels configured  (four on my machine and four on a machine in computer lab) plus a master or controlling kernel.  The controlling kernel has a KernelID of zero and these have ID's of one through eight. I can check this by evaluating $KernelID on each of the kernels

ParallelEvaluate[$KernelID]

Or by looking at the parallel kernel status window which you will find under the Evaluation menu. Mine is shown in the figure below.

The communication overhead can be checked by asking a kernel for it's KernelID.  Checking a local kernel on my machine

t1=AbsoluteTime[];
ParallelEvaluate[$KernelID,6];
timeLocal=AbsoluteTime[]-t1

and comparing the response time with a kernel on the machine in the computer lab

t1=AbsoluteTime[];
ParallelEvaluate[$KernelID,1];
timeLab=AbsoluteTime[]-t1

shows that it takes almost 50% longer. The machine in the lab is on a different sub-net than mine but the network between them is fast with an average ping time to the machine (at the time of the test) of less then 1/2 a millisecond.

timeLab/timeLocal

t1=AbsoluteTime[];
$KernelID;
timeMaster=AbsoluteTime[]-t1

and both take longer than the same request on the master kernel.

timeLocal/timeMaster

Labs

Computer labs are easy to turn into a grid using Mathematica's Lightweight Grid Services. The disadvantage of using a computer lab as grid resource is that communication between machines is relatively slow as we saw in the example above since the machines are just connected by standard Ethernet. The machines may also not all be the same leading to an inhomogeneous grid. Depending on the type of calculation you are doing inhomogeneous grids may or may not be a problem.

Dedicated Clusters

Dedicated cluster range from small Beowulf clusters to the clusters at national supercomputer centers. Big clusters offer a host of advantages, they may have thousands, tens of thousands or even a two hundred thousand cores (Blue Waters), very fast interconnects between the cores, abundant memory and large fast file systems.
The disadvantage is that they are harder to use. You must apply for an allocation and learn to use the clusters queue system. Jobs of any significance need to be submitted as batch jobs and this makes debugging much harder. You should always try to debug on a small local grid where you can run interactively before submitting a batch job.

Code Performance

Predicting how much speed-up you will get by running on multiple cores and how this will scale with the number of cores is very difficult. It depends not only on the percent of your code that can be run in parallel and on the latency between processors but also on factors such a cache size, chip architecture and of course your algorithm. A better algorithm will almost always beat more cores! It is always a good idea to run some benchmarks for your problem before embarking on a large run.

Parallel Computation

The details of parallel computation depend of course on the problem but roughly speaking problems fall into two classes, those that are embarrassingly parallel and those that are not. Embarrassingly parallel problems (now sometimes called perfectly parallel problems for PR reasons) are those where every processor does an independent calculation. Some examples will make the distinction clear.

Examples

Table or DO loop in Parallel

As simple and common case is where we have a nested Table command or Do loop which we want to split up over a number of processors. The command for this in Mathematica is ParallelTable. Let's look at some examples. Suppose we want to compute 16 million random numbers. We can do this a single processor with

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

and it only takes 1.7 seconds.

Seeking to speed things up I launch four kernels on my machine and use ParallelTable.

ParallelEvaluate[$KernelID]

resulting in a slow down by a factor of more than 2.6.

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

For parallel computation to be useful the computation time must dominate the time spent by the processors communicating with one other. In this example the computation is trivial but each processor has to send two million numbers to the controlling processor. So let's try a case that's more computational intensive. In this example we evaluate a set of integrals that arise the calculation of the of bound states for a particle partially confined to an irregularly shaped region.  Each integral takes a second or so to compute and an accurate calculation of the bound states requires over 20,000 such integrals. That makes this a ideal problem for grid computing, it's computational intensive, the calculations are independent and very little information needs to be exchanged between the processors. Run on just on processor our sample calculation takes 283 seconds. You will not get exactly the same result each time because the integration method used is a Monte Carlo method.

t1=AbsoluteTime[];results=Table[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"],{i,1,8},{j,1,8}]
timing1=AbsoluteTime[]-t1

However if we run on eight processors (four local and four on a machine in our computer lab it only takes 50 seconds, a speed-up of more than 5 and a half times.

t8=AbsoluteTime[];results=ParallelTable[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"],{i,1,8},{j,1,8}]
timing8=AbsoluteTime[]-t2

Shown below is a picture the showing the activity on the parallel kernels.

Now suppose we want to run on 16 processors. I have added a couple of additional machines to my cluster as show by the list of kernel id's.

ParallelEvaluate[$KernelID]

We cannot simply use ParallelTable anymore since that will distribute the computations over only as many processors as the length of the outer list which is this case is eight. Instead I will generate a list of all the indices that I want to evaluate the integral over

data=Table[{i,j},{i,1,8},{j,1,8}]

and define a function that takes the an ordered pair and returns the numerical integral associated with that order pair.

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"]

The definition of the function needs to be known on all the kernels and this is done using the command DistributeDefinitions.

DistributeDefinitions[f]

We can now compose the function f with ParallelSubmit using the function Composition. This can then be mapped (see Map) across the flattened list of indices. This sends each computation off to the next available kernel. The function WaitAll waits until all of the computations are finished and then returns the result.

t16=AbsoluteTime[];result=WaitAll[Map[Composition[ParallelSubmit,f],Flatten[data,1]]]
AbsoluteTime[]-t16

You will notice that we get a speed-up of almost but not quite a factor of two over using eight processors. If we like we can recover the array by partitioning the list into groups of eight.

Partition[result,8]

To get a better feeling for the scaling I will add a few more machines to my cluster so that I have 32 processors

ParallelEvaluate[$KernelID]

and re-run

t32=AbsoluteTime[];result=WaitAll[Map[Composition[ParallelSubmit,f],Flatten[data,1]]]
AbsoluteTime[]-t32

Plotting the speedup versus the number of processors shows that although we get a substantial speedup the scaling is not linear. That is using n processors does not produce a speedup by a factor of n. Since there is always some communication overhead you can never do as well as linear scaling.

Monte Carlo in Parallel

Monte Carlo calculations form an important class of problems that are embarrassingly parallel. Suppose we want to calculate π which is a classic Monte Carlo problem. We can calculate π by generating a list of ordered pairs of random numbers that lie in the unit rectangle. We then check to see which of these lie with in the unit circle and divide the number in the unit circle by the total number. This will give π/4. The geometry is shown below for a random set of 200 ordered pairs.

points=Table[{Random[]^2,Random[]^2},{i,1,200}];

Show[Graphics[{Red,Circle[{0,0},1]}],Graphics[{Blue,Opacity[.5],Rectangle[{0,0},{1,1}]}],Graphics[Map[Point,points]]]

To check to see if the ordered pair lies inside the unit circle we compute the distance d the point is from the unit circle and them compute 1-Floor[d]. Since 0≤ d≤ (Sqrt[2]), 1-Floor[d] will give zero if d>1 and 1 if d<1. Summing the result then gives the number in the that lie inside the unit circle. The code below computes π on one processor using  40,000 points.

t1=SessionTime[];
m=10000;
n=4;
acceptpoint[j_]:=Total[Table[1-Floor[Sqrt[(Random[])^2+(Random[])^2]],{i,1,m}]]
hits=Table[acceptpoint[j],{j,1,n}];
hits=Total[hits];
pi=hits/(n*m)*4.0
t2=SessionTime[]-t1

We now repeat the calculation with the number of points running from 40000 to 2400000.

data=Table[t1=SessionTime[];
n=4;
m=10000 i;
acceptpoint[j_]:=Total[Table[1-Floor[Sqrt[(Random[])^2+(Random[])^2]],{i,1,m}]];
hits=Table[acceptpoint[j],{j,1,n}];
hits=Total[hits];
pi=hits/(n*m)*4;
t2=SessionTime[]-t1,{i,1,60,1}]

size=Table[10000 n i,{i,1,60,1}];

and plot the timings.

plot1=ListPlot[Transpose[{size,data}],PlotStyle->{Red,PointSize[0.02]},AxesLabel->{"number of trials","Wall Time"}]

I now launch four local kernels on my machine  and run the in parallel by using ParallelTable instead of Table.

ParallelEvaluate[$KernelID]

paralleldata=Table[t1=SessionTime[];
n=4;
m=10000 i;
acceptpoint[j_]:=Total[Table[1-Floor[Sqrt[(Random[])^2+(Random[])^2]],{i,1,m}]];
DistributeDefinitions[n,m,acceptpoint];
hits=ParallelTable[acceptpoint[j],{j,1,n}];
hits=Total[hits];
pi=hits/(n*m)*4;
t2=SessionTime[]-t1,{i,1,60,1}]

and plot the timings.

plot2=ListPlot[Transpose[{size,paralleldata}],PlotStyle->{Blue,PointSize[0.02]},AxesLabel->{"number of trials","Wall Time"}]

Combining the plots shows that we get a significant speed up. The difference between this and the early table example is that here each processor sends only a single number to the master kernel, not a large table so there is very little communication overhead.

Show[plot1,plot2,ImageSize->72 6,PlotLabel->"Red = Single Processor \n Blue = Four Processors"]

Simulated Annealing

Although it is based on Monte Carlo methods simulated annealing is not embarrassingly parallel since the state of the whole system is needed at each step. There are a number of algorithms for parallel simulated annealing and choosing the best one for your problem in non-trivial[2,3].

Algorithm testing

Another application of parallel computing is algorithm testing. A classic case is global minimization where there is no method that is guaranteed to find the global minimum in a finite amount of time. In Mathematica global minimization is done by the function NMinimize. NMinimize can use four different method, differential evolution, Nelder-Mead, random search and simulated annealing. It is often very difficult to know for a given problem which method is best. With parallel computing we can test them all at once. 

A test problem given in the Mathematica documentation is to find the global minimum of the function f=(1/4) ((x^2)+(y^2))-sin(10 (x+y))+(e^(sin(50 x)))+sin(70 sin(x))+sin(60 (e^y))+sin(sin(80 y)). Plotting shows that the function has a huge number of local minima.

Clear[f];
f=(E^(Sin[50 x]))+Sin[60 (E^y)]+Sin[70 Sin[x]]+Sin[Sin[80 y]]-Sin[10 (x+y)]+(1/4) ((x^2)+(y^2));
Plot3D[f, {x, -1, 1}, {y, -1, 1}, PlotPoints->50, Mesh->False]

There are four possible methods so I will launch four kernels

ParallelEvaluate[$KernelID]

and define a list with the four methods as strings.

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

I need to distribute the definition of the function f.

DistributeDefinitions[f]

I can now use ParallelMap try each method on a different kernel.

ParallelMap[{$KernelID,#,Timing[NMinimize[f,{x,y},Method->#]]}&,methods]

We see that differential evolution, run on kernel number 4, is by far the slowest but returns the best answer. If we let NMinimize use the default method it uses Nelder-Mead which is fast but gives only the third best answer on this problem.

NMinimize[f,{x,y}]

The Diffusion Equation

The numerical solution of partial differential equations is one of the most important topics in computational physics. As an example consider the heat or diffusion equation ((∂U)/(∂t))=κ(((∂^2)U)/(∂(x^2) )). One way to solve this equation is to put the spatial part on a grid and then step forward in time.  This easy to implement but is only conditional stable, we must take time steps no larger than ((Δx)^2)/(2κ) where Δxis the spatial step in order to avoid instabilities. A better method which is unconditionally stable is the Crank-Nicolson method [4]. In the Crank-Nicolson method the spatial derivatives are discretized using a central difference and the time derivative is approximated with the trapezoidal rule. That means that ((∂f(t, x))/(∂t))=g(t,x)->((f((Subscript[t, j+1]), (Subscript[x, i]))-f( (Subscript[t, j]), (Subscript[x, i])))/(Δt))=(1/2)(g((Subscript[t, j+1]),(Subscript[x, i]))-g((Subscript[t, j]),(Subscript[x, i]))). The advantage of this method is that it is unconditionally stable so we can take much larger time steps. The disadvantage is that the (Subscript[t, j+1])terms appear on both sides of the equation so we must solve a set of linear equations at each time step. Using the trapezoidal rule for the time derivative and central differences for the spatial derivatives converts the diffusion equation into 
 ((u((Subscript[t, j+1]), (Subscript[x, i]))-u((Subscript[t, j]), (Subscript[x, i])))/(Δt))=(κ/(2((( Δx))^2)))(u((Subscript[t, j+1]),(Subscript[x, i+1]))-2u((Subscript[t, j+1]),(Subscript[x, i]))+u((Subscript[t, j+1]),(Subscript[x, i-1]))+u((Subscript[t, j]),(Subscript[x, i+1]))-2u((Subscript[t, j]),(Subscript[x, i]))+u((Subscript[t, j]),(Subscript[x, i-1]))).
 It is customary to define r=((κ Δt)/(((Δx))^2))giving -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])).
If we take r=1 the equations simplify to
 -u((Subscript[t, j+1]),(Subscript[x, i+1]))+4u((Subscript[t, j+1]),(Subscript[x, i]))- u( (Subscript[t, j+1]),(Subscript[x, i-1]))  =
u( (Subscript[t, j]),(Subscript[x, i+1]))+ u((Subscript[t, j]),(Subscript[x, i-1])). 
We could use larger values for r, this would allow larger time steps. There would be some loss of accuracy but the method would still be stable.

To implement the Crank-Nicolson method we start by defining the initial concentration. To begin with I will work on a 2000 point grid. So if the grid covers say a meter then Δx=0.0005 meters and thus κ Δt=2.5 (10^(-7))(m^2)/s. Thus if κ=1(m^2)/s each time step is 2.5 (10^(-7))s which is very small. We could get larger time steps by using a larger value of r.

InitialConcentration=Table[((N[Sin[2Pi i/2000]])^2),{i,0,2000}];
plotInitial=ListPlot[InitialConcentration,Joined->True]

We need to impose the boundary conditions which we will take to be sinks and initialize u(0,x).  We also need to define the variables u(j,i) that we will solve for at each time step j.

Clear[u,r]
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]}]

We can now define the equations for Crank-Nicholson which we will do for arbitrary r. The function "equations[j]" generates the full set of equations for time step j.

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]}]

The equations that we need to solve are of the form M x=b where M is tridiagonal. Taking advantage of Mathematica's SparseArray and PackedArray technologies we can quickly solve very large systems of tridiagonal equations. We first use CoefficientArrays to generate the matrix M and vector b. This can be slow for large arrays but we only need to do it once since M is the same at each time step. We will however have to recompute b with each time step.

r=1;
t1=TimeUsed[];
{b,M}=-CoefficientArrays[equations[0],var[1]]
TimeUsed[]-t1

Notice that both M and b have been returned as SparseArray objects which take very little storage. We can check this by using ByteCount

ByteCount[M]

which shows that M takes up less than 81 kilobytes but the normal form of M takes just over 96 megabytes.

ByteCount[Normal[M]]

ByteCount[M]/ByteCount[Normal[M]]//N

To compute b at each step j we need u[j, i - 1] + u[j, i + 1] + (2/r - 2) u[j, i]. We can efficiently generate the u[j, i - 1] + u[j, i + 1]  terms by list rotation (see RotateRight and RotateLeft) as the small symbolic example below shows.

testData=Table[U[i],{i,1,5}]

RotateRight[testData]+RotateLeft[testData]

This gets the interior elements right but not the ends. To fix the ends we insert (see Insert) the left and right boundary conditions (note the factors of two) at the start and end of the list

RotateRight[Insert[Insert[testData,2leftend,1], 2rightend,-1]]+RotateLeft[Insert[Insert[testData, 2leftend,1],2 rightend,-1]]

and then drop the first and last elements

Take[RotateRight[Insert[Insert[testData,2leftend,1], 2rightend,-1]]+RotateLeft[Insert[Insert[testData, 2leftend,1],2 rightend,-1]],{2,-2}]

The  (2/r - 2) u[j, i] term is then just added to this result.  Putting this together gives

Do[newU=LinearSolve[-M,b];
b=Take[RotateRight[Insert[Insert[newU,2leftend,1], 2rightend,-1]]+RotateLeft[Insert[Insert[newU, 2leftend,1],2 rightend,-1]],{2,-2}]+(2/r-2) newU,{i,1,6000}]//Timing

Since LinearSolve uses very efficient methods for tridiagonal systems it takes only 5.7 seconds to do 6000 iterations.

ListPlot[{InitialConcentration,newU},PlotLabel->"Initial concentration  in blue \n concentration after 6000 steps in red",ImageSize->6 72 ]

If we change r we take large time steps (although we lose some accuracy). If for example we take r=10 we need only 600 steps.

r=10;
{b,M}=-CoefficientArrays[equations[0],var[1]];

Do[newULargeR=LinearSolve[-M,b];
b=Take[RotateRight[Insert[Insert[newULargeR,2leftend,1], 2rightend,-1]]+RotateLeft[Insert[Insert[newULargeR, 2leftend,1],2 rightend,-1]],{2,-2}]+(2/r-2) newULargeR,{i,1,600}]//Timing

As you can see the change is small.

ListPlot[newU-newULargeR,PlotRange->All]

To solve the problem in parallel we break the spatial domain into pieces and give each processor one or more pieces. This problem is not embarrassingly parallel since at the end of each time step the processors have to exchange information about the boundaries of their pieces.

Information exchanged between processors in domain decomposition.

There are two general approaches for this. Synchronous methods in which each processor communicates the current values of it's boundaries to it's neighboring processors and asynchronous methods in which the boundary variables are allowed to lag by one time step. Synchronous methods are more accurate but slower. Asynchronous methods are faster and easier to program but make (hopefully small) errors at the boundaries. We will  look at an asynchronous method.

An Asynchronous method

In order to develop the code we start with a small grid of forty points.

InitialConcentration=Table[((N[Sin[2Pi i/40]])^2),{i,0,40}];
plotInitial=ListPlot[InitialConcentration,Joined->True]

The initialization as the same as for a single processor.

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]}]

As is the equation generation.

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]}]

We first set the number of processors (in this example four) and break the list into the same number of parts so that each part can be sent to a processor.

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]}]

We next map CoefficientArrays over the segments to generate the matrices and vectors that we will send to LinearSolve.

r=10;arrays=Map[CoefficientArrays[#1[[1]],#[[2]]]&,segments[1]];
vectors=-Map[First,arrays];
matrices=Map[Last,arrays];

Checking shows that we do indeed have four different matrices and four different vectors

matrices

vectors

All the processors need to know the definitions of the matrices and vectors so we need to use DistributeDefinitions to achieve that.

DistributeDefinitions[matrices,vectors]

Each one of our vectors contains an undefined variable. In the first vector this is u(1,12) which lies on the left edge of the second vector. You can check this by converting the first vector from a sparse array to a normal object.

vectors[[1]]//Normal

The first line of code in the next cell u[1,x_]=u[0,x] will set u[1,12] in vector one to u[0,12]. This is why the code is called asynchronous.  This leads to an error but the error will be small if the time step is sufficiently small. The next line of code ships LinearSolve[matrices[[i]],vectors[[i]]] off to processor i (see MapThread).

u[1,x_]=u[0,x];
newU=WaitAll[MapThread[ParallelSubmit[LinearSolve[#1,#2]]&,{matrices,vectors}]];

Next, we must find the values of the vectors at the edges of each domain

edgeValues=Map[Drop[#,{2,-2}]&,newU];

edgeValues

and impose the boundary conditions and the left and right most edges (see ReplacePart)

boundaries=ReplacePart[edgeValues,{{1,1}->leftend,{n,2}->rightend}]

We can bundle this all up in Do loop. The only new wrinkle is that we have to use MapThread to form the new vectors since we now have as many vectors as processors. We can test the speed and accuracy of our code with just a few iterations (in this case four)

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,4}]
t2=AbsoluteTime[]-t1

In order to plot the new concentration we need to flatten the list since the result is returned as a nested list with as many sublists as processors.

plotMultiple=ListPlot[Flatten[newU],Joined->True]

We need to check our result against the single processor result and the code for a single processor is bundled in the cell below.

Clear[u,b,M]
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]}]

{b,M}=-CoefficientArrays[equations[0],var[1]];

t1=TimeUsed[];
Do[newUSingle=LinearSolve[-M,b];
b=Take[RotateRight[Insert[Insert[newUSingle,2leftend,1], 2rightend,-1]]+RotateLeft[Insert[Insert[newUSingle, 2leftend,1],2 rightend,-1]],{2,-2}]+(2/r-2)newUSingle,{i,1,4}];
TimeUsed[]-t1

plotSingle=ListPlot[newUSingle]

Comparing the results shows that our multiple processor code took 25 times as long to run and returned an inaccurate answer. But all is not lost. The problem is our small grid. As small grid is good for debugging but you would never run a problem this small on multiple processors. The task sent to each processor, solving 10 linear equations is trivial so the problem is dominated by the communications overhead. Furthermore we are taking very large time steps. Since  r=((κ Δt)/(((Δx))^2))and we took r=10 we have Δt=10 (((Δx))^2)/κ. If κ=1(mm^2)/s and our spatial domain is 1cm then with forty grid points Δx= 0.25mm and Δt=0.625 seconds.

Show[plotMultiple,plotSingle]

If we try a larger grid things look better. If we take a grid of 40000points

InitialConcentration=Table[((N[Sin[2Pi i/40000]])^2),{i,0,40000}];
plotInitial=ListPlot[InitialConcentration,Joined->True]

and run for 4000 steps (the code for this is given in the next cell)

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

and compare this to the single processor case

Clear[u,b,M]
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]}]

{b,M}=-CoefficientArrays[equations[0],var[1]];
t1=TimeUsed[];
Do[newUSingle=LinearSolve[-M,b];
b=Take[RotateRight[Insert[Insert[newUSingle,2leftend,1], 2rightend,-1]]+RotateLeft[Insert[Insert[newUSingle, 2leftend,1],2 rightend,-1]],{2,-2}]+(2/r-2)newUSingle,{i,1,4000}];
TimeUsed[]-t1

we see that we are now only 1.5 times slower (not 25) and accuracy is good. You can clearly see the errors (which are now small) at the edges of the domains due to the asynchronous updating.

ListPlot[Flatten[newU]-newUSingle,PlotRange->All,PlotLabel->"Error due to  asynchronous updating",AxesLabel->{"grid point", "Difference"}]

The plot below shows the speed up or slow down as a function of the grid size.

How much speed up you can expect as you increase the number of processors? That depends on the size of your problem and the architecture of your grid.  The plot below shows timings for a grid of one million elements run on Glenn at OSC. The code was run on dual socket quad core Opterons. Notice that using six processors is  slower than using four that because the code is now split over two chips not just cores on the same chip so the  I/O is slower..

Running on a large machine

As an example I have reproduced the files I used to run on the Ohio Supercomputer Center's (OSC) supercomputer Glenn. On a large cluster you will not be running interactively but will have to submit a batch job. The details will depend on the batch system used on the machine you are running on. OSC uses PBS. Information on using PBS (Portable Batch System) can be found here.

Set Up

The first step is to create a batch file. The first line tells the scheduler that I want up to 50 hours of wall time (that is the time shown by a clock on the wall). The next line asks for five nodes with four processors per node. The third line says that my job will be called nef (for numerical eigenfunctions). The line #PBS -j oe says that output and error messages should be written to standard out and standard error. The #PBS -S/bin/ksh makes sure that I am in the Korn shell. The cd $Home command changes the directory to my home directory just in case I am not already there. The last line gives the file path tho Mathematica 6.02 on Glenn (this will be different on other machines),  launches the kernel without a prompt (you don't want a prompt in a batch job), and reads in the file NEF.m.

Create a file called my_request.job with the following code

# PBS - l walltime = 50 : 00 : 00
# PBS - l nodes = 5 : ppn = 4
# PBS - N nef
# PBS - j oe
# PBS - S/bin/ksh
cd $Home
/usr/local/Wolfram/Mathematica7 .0 .1/OSC/math - noprompt - run "<<NEF.m"

The file NEF.m

The file is reproduced below. Those lines highlighted in yellow are the changes I have made to run the code in parallel.  The first line loads the cluster integration package. The next line tells the cluster engine that I will be using PBS with the torque scheduler. The option MathKernelCommand->env OMP_NUM_THREADS=1 tells the Intel math libraries not to thread processes over multiple cores (we will be doing this ourselves). Setting the ToQueue option to True simply says that I will be using the queue. I don't need a queue name since I have a batch file that handles the naming for me. After the cluster setup I need to add one line of code and change one line. The line added is

DistributeDefinitions[L, hbar, HK, phi, integral, m, Voutside, Vdiff, top, bottom, xedge, yedge, topboundary, bottomboundary, a, numberOfBasisFunctions];

which exports the definitions to all the worker kernels. You can read about resource sharing in parallel computing here. The other change is to change Table to ParallelTable in the computation of the numerical integrals. I have 18 basis function so I need to do (18^4)=104,976  numerical integrals. Since I have 18 processors I use parallel table to ship each processor (18^3) integrals. At the end of the computation I write the results to a file which I tag with the date and time so that I can keep the files straight.

Needs["ClusterIntegration`"];
ClusterSetEngine[PBS -> {
        "ClusterEnginePath" -> "/usr/local/torque-2.1.8", "MathKernelCommand" -> "env OMP_NUM_THREADS=1 /usr/local/Wolfram/Mathematica7.0.1/OSC/math",
        "ToQueue" -> True,
        "NativeSpecification" -> {"QueueName" -> None}}
    ];
CPUs = 18;
ClusterLaunchSlaves[CPUs];
Clear[L];{xedge,yedge}={x,y}/.NMaximize[{y,x^2+Sin[x y]+ y^2/4==1},{x,y}][[2]];
topboundary[y_/;0<=y<=yedge]:=x/.FindRoot[x^2+Sin[x y]+ y^2/4==1,{x,0,xedge,1.5}];
topboundary[y_/;y<0]:=-bottomboundary[-y];
bottomboundary[y_/;0<=y<=yedge]:=x/.FindRoot[x^2+Sin[x y]+ y^2/4==1,{x,-1.,-1.5,xedge}];
bottomboundary[y_/;y<0]:=-topboundary[-y];
Quiet[
top=FunctionInterpolation[topboundary[y],{y,-yedge,yedge}];
bottom=FunctionInterpolation[bottomboundary[y],{y,-yedge,yedge}]
];
Voutside[x_,y_]=(1/2)(x^2+y^2)+100;
Vdiff[x_,y_]=-100;
integral[phi1_,phi2_,options___]:=NIntegrate[phi1 Voutside[x,y]phi2,{y,-L/2,L/2},{x,-L/2,L/2},options]+NIntegrate[phi1 Vdiff[x,y]phi2,{y,-yedge,yedge},{x,bottom[y],top[y]},options];
timeone=AbsoluteTime[];
phi[i_,j_,x_,y_]=((Sqrt[2/L] Sin[i Pi ( x+L/2)/L])(Sqrt[2/L] Sin[j Pi ( y+L/2)/L]));
HK[i_,j_,k_,l_]=Integrate[(phi[i,j,x,y])(-hbar^2/(2m)( D[phi[k,l,x,y],{x,2}]+D[phi[k,l,x,y],{y,2}])),{x,-L/2,L/2},{y,-L/2,L/2},Assumptions->{Element[i,Integers],Element[j,Integers],Element[k,Integers],Element[l,Integers]}];
HK[i_,j_,i_,j_]=Integrate[(phi[i,j,x,y])(-hbar^2/(2m)( D[phi[i,j,x,y],{x,2}]+D[phi[i,j,x,y],{y,2}])),{x,-L/2,L/2},{y,-L/2,L/2},Assumptions->{Element[i,Integers],Element[j,Integers],L>0}];
L=10;
hbar=1;
m=1;
numberOfBasisFunctions=CPUs;
DistributeDefinitions[L,hbar,HK,phi,integral,m,Voutside,Vdiff,top,bottom,xedge,yedge,topboundary,bottomboundary,a,numberOfBasisFunctions];
a[i_,j_,k_,n_]:=HK[i,j,k,n]+integral[phi[i,j,x,y],phi[k,n,x,y],AccuracyGoal->3,PrecisionGoal->3];
M=ParallelTable[a[i,j,l,n],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions},{l,1,numberOfBasisFunctions},{n,1,numberOfBasisFunctions}];
M1=Partition[Flatten[M],numberOfBasisFunctions^2];
results=Chop[Eigenvalues[M1], 10^-9];
basisFunctions=Flatten[Table[phi[i,j,x,y],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions}]];
wavefunctions=Eigenvectors[M1].basisFunctions/Norm[Eigenvectors[M1].basisFunctions];
Put[results,StringJoin["energies",ToString[Date[][[3]]],"_",ToString[Date[][[4]]],"_",ToString[Date[][[5]]]]];
wavefunctionData=Take[wavefunctions,{Length[wavefunctions]-9,Length[wavefunctions]}];
Put[wavefunctionData,StringJoin["wavefunctions",ToString[Date[][[3]]],"_",ToString[Date[][[4]]],"_",ToString[Date[][[5]]]]];
Quit[];

Now of course we need to submit the job. This is done with the command (from a UNIX shell logged into the cluster) 
qsub my_request.job
qsub stands for queue submit.

If you try this project it is critical that you first develop code that works on a single processor for a small number of basis functions and then port the code to a big machine!

References

[1]	Rubin Landau, Manuel José Páez and Cristian C. Bordeianu, A Survey of Computational Physics, page 36,  Princeton University Press, 2008
[2]	DJ Chen, CY Lee, CH Park, P Mendes, Parallelizing simulated annealing algorithms based on high-performance computer, Journal of Global Optimization, 39, 261-289,  2007 - Springer
[3]	Daniel R Greening, Parallel Simulated Annealing Techniques, Physica D 42, 293-3006, 1990
[4]	J. Crank and P. Nicolson, A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type, Proc. Cambridge Phil. Soc 43, 50, 1946.