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)