# Lock-free cache oblivious n-body algorithm

I'm currently looking at, from a rather high level, the parallelization of the gravity calculation in an N-body simulation for approximating a solution to the N-body problem.

The simple form of the algorithm looks something like this:

for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
if (i != j)
bodies[i].ApplyGravityFrom(bodies[j]);

Body::ApplyGravityFrom(Body &other)
{
Vector dr = other.Pos - this->Pos;
double r2 = Dot(dr, dr);
double ir3 = 1 / (r2 * sqrt(r2));
this->Acc += (other.Mass * ir3) * dr;
}


This simple version has an obvious parallelization over the outer loop:

#pragma omp parallel for
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
if (i != j)
bodies[i].ApplyGravityFrom(bodies[j]);


However, you're doing twice as many units of work as necessary. It's the case that gravity acting on body i from body j is the same as the gravity acting on body j from body i, but with the opposite sign.

You can calculate gravity pairwise instead:

for (int i = 0; i < n; i++)
for (int j = i + 1; j < n; j++)
bodies[i].PairwiseGravity(bodies[j]);

Body::PairwiseGravityBody &other)
{
Vector dr = other.Pos - this->Pos;
double r2 = Dot(dr, dr);
double ir3 = 1 / (r2 * sqrt(r2));
this->Acc += (other.Mass * ir3) * dr;
other.Acc -= (this->Mass * ir3) * dr;
}


This is the same exact calculation but you're making use of the fact that the force is symmetric, but with a sign flip.

But now parallelizing this loop becomes much harder. But if you treat this problem geometrically, then you can realize that you can actually divide up the work in a special way to never require locks.

First off, you can calculate pairwise forces for the first n / 2 bodies. This is inherently serial. Then, you can calculate the pairwise forces between the bodies [0, n / 2] and [n / 2, n]. Lastly you calculate the pairwise forces for the last n / 2 bodies.

The algorithm then looks like:

void BlockGravity(Body *bodies, int i1, int i2, int j1, int j2)
{
for (int i = i1; i < i2; i++)
for (int j = j1; j < j2; j++)
bodies[i].Pairwise(bodies[j]);
}

void TriangleGravity(Body *bodies, int i1, int i2)
{
for (int i = i1; i < i2; i++)
for (int j = i1 + 1; j < i2; j++)
bodies[i].Pairwise(bodies[j]);
}

void Gravity(Body *bodies, int n)
{
TriangleGravity(bodies, 0, n / 2);
BlockGravity(bodies, n / 2, n, 0, n / 2);
TriangleGravity(bodies, n / 2, n);
}


The block gravity calculations have the property that none of them ever write into the same memory twice. Every pair used in that calculation is unique. What's special is at this point you can actually continue to subdivide the problem further:

void Gravity(Body *bodies, int n1, int n2)
{
int n = n2 - n1;
if (n > MinBlockSize)
{
int m = n1 + n / 2;
Gravity(bodies, n1, m);
Gravity(bodies, m, n2);
BlockGravity(bodies, m, n2, n1, m);
}

else
TriangleGravity(bodies, n1, n2);
}


By subdividing the problem, you gain the benefit of cache locality. The smaller problems tend to work on things that close to each other, whereas the naive version would continuously flush the cache as it iterated over the whole array. I included a minimum size for subdivision since at some point subdividing the problem any further won't get any benefits

So far, this does the same amount of work as simply calling TriangleGravity(bodies, 0, n). There is an opportunity for parallelism in the recursive calls to Gravity, as well as the block computation. I'll post this once I fully test it.

Does there seem to be any logical flaws in this? I'll try to throw up an image visually demonstrating what it does, because seeing it on paper helps A LOT.

Update

Here's some performance data for the scaling with adjusting the MinBlockSize parameter. These tests were performed for n = 4096 with fully parallelism being exploited within the recursive calls.

BlockSize \ Timing (milliseconds)
1    - 201 (Overkill)
2    - 144
4    - 123
8    - 106
16   - 101
32   - 110
64   - 101 (Optimal Performance)
128  - 108
256  - 124
512  - 147
1024 - 195
2048 - 293
4096 - 407 (Same as naive)


You seem to be making the problem more complex:

I would take a step back to your original algorithm:

#pragma omp parallel for
for (int i = 0; i < n; i++)
for (int j = i + 1; j < n; j++)
bodies[i].PairwiseGravity(bodies[j]);


Your problem (as you stated) that the parallelism is uneven (the lower outer loops is doing more work than the higher values of the outer loop).

But why not combine the two loops.

#pragma omp parallel for
for (int i = 0; i < ((n*(n-1))/2; i++)
bodies[src(i,n)].PairwiseGravity(bodies[dst(i,n)]);


Each loop is completely independent.
All you have to do is calculate how the functions src() and dst() work.

• The code I posted improves performance significantly without parallelism. For n = 1024, the naive algorithm takes approximately 2.5 milliseconds to execute. Using MinBlockSize = 64 the same calculation takes 0.84 milliseconds to execute. The cache effects aren't to be underestimated. I added the definition of PairwiseGravity, and you'll see the only issue with what you suggest is that it would be extremely easy to have two threads writing into the same memory at once. – Mike Bailey Jan 8 '12 at 17:03
• If for n = 1024 and time of 2.5 milliseconds you probably should not even be considering threading. But your technique stills seems a bit costly (in terms of maintenance and initial creation) for such a small gain (admittedly 3 fold increase is good but will it still scale for a data set size that would actually take significant time). If you actually did some real work (meant in terms of the functions work(Not a comment on your code)) then you should be considering a technique called map reduce. – Martin York Jan 9 '12 at 3:12
• That's a tiny problem size. Most things I'm working on are for far larger n, and this grows to O(n^2) so a more realistic n = 10000 would take 250 milliseconds per iteration. On my application currently, I've actually managed to get fully linear speed up with every core added. A previous version with a simpler threading scheme had worse performance and started degrading as the problem size increased. One of the big issues is the data stored in these simulations isn't small -- it's on the order of 1 KB per body (e.g. n = 10,000 -> 10 MB data) – Mike Bailey Jan 9 '12 at 3:14
• My main reason for posting here was to figure out if I did this appropriately :) I actually had a bug in my implementation that was posted, I'll update this shortly to account for that. – Mike Bailey Jan 9 '12 at 3:16
• So you are threading I thought you just said you were not! OK. Now that we are getting down to it. Maybe if you provided more detail about the actual work being done then we can consider a proper technique like Map Reduce. – Martin York Jan 9 '12 at 3:19

There is no logical flow; this approach is known, and indeed it benefits from better cache locality, and is very suitable for parallelism, especially with work-stealing-based frameworks.

Note that you can calculate forces between the first N/2 bodies and the last N/2 bodies (i.e. the triangles) in parallel, then you might process the block. When processing the block, you might also apply recursion and split it in 4 sub-blocks (rather than by rows or cols). These 4 sub-blocks can be split into two pairs along diagonals of the original block; in each pair, sub-blocks don't share data and can be processed in parallel:

void BlockGravity(Body *bodies, int i1, int i2, int j1, int j2)
{
int i_mid=(i1+i2)/2, j_mid=(j1+j2)/2;
{ // This pair can be processed in parallel
BlockGravity(bodies, i1, i_mid, j1, j_mid);
BlockGravity(bodies, i_mid, i2, j_mid, j2);
}
// Synchronize here
{ // This pair can also be processed in parallel
BlockGravity(bodies, i1, i_mid, j_mid, j2);
BlockGravity(bodies, i_mid, i2, j1, j_mid);
}
}


Another description of this algorithm (with pictures) is here: http://software.intel.com/en-us/blogs/2010/07/01/n-bodies-a-parallel-tbb-solution-parallel-code-a-fresh-look-using-recursive-parallelism/, and subsequent posts in that series contain a parallel implementation sketch using Intel's TBB as well as performance measurements.