Lock-free cache oblivious n-body algorithm

Question

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)

Answer 1

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.

Answer 2

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.