The n-Body Problem and Shared-Memory Solvers

In the preceding chapters, we explored fundamental parallel programming APIs in isolation: MPI for distributed-memory systems, OpenMP and Pthreads for shared-memory architectures, and CUDA for massively parallel GPUs. In this final chapter, we transition from individual APIs to end-to-end parallel program development.

We will examine two realistic computational problems:

1. The $n$-body problem: Simulating the gravitational interactions of celestial bodies or molecular particles over time.
2. Sample sort: A scalable parallel sorting algorithm designed for unknown data distributions.

By designing solutions across all three architectural paradigms, we will uncover striking similarities in work partitioning, identify patterns where parallel algorithms diverge sharply from their serial counterparts, and learn how to select the right parallel API for a given problem.

---

7.1 The $n$-Body Problem

In an $n$-body problem, we seek the positions and velocities of a collection of $n$ interacting particles over a specified duration of time. Applications span astrophysics (simulating stars, galaxies, and planetary systems) and molecular dynamics (modeling protein folding and atomic interactions).

The input to the simulation comprises:

• The number of particles $n$.
• The mass $m_q$, initial position $\mathbf{s}_q(0)$, and initial velocity $\mathbf{v}_q(0)$ for each particle $q \in \{0, 1, \dots, n-1\}$.
• The timestep length $\Delta t$ and the total simulation duration $T \Delta t$.

The output consists of the positions and velocities of all $n$ particles at discrete timesteps $t = 0, \Delta t, 2\Delta t, \dots, T\Delta t$, or simply at the final time $T\Delta t$. To maintain computational clarity while preserving the mathematical essence, we assume motion occurs in a two-dimensional plane ($\text{DIM} = 2$).

---

7.1.1 Physical Laws and Mathematical Formulation

The physical model is governed by Newton’s second law of motion ($\mathbf{F} = m \mathbf{a}$) and his law of universal gravitation.

The gravitational force $\mathbf{f}_{qk}(t)$ exerted on particle $q$ by particle $k$ at time $t$ is:

$\mathbf{f}_{qk}(t) = -\dfrac{G m_q m_k}{|\mathbf{s}_q(t) - \mathbf{s}_k(t)|^3} \left[ \mathbf{s}_q(t) - \mathbf{s}_k(t) \right]$

Where:

• $G$ is the gravitational constant ($6.673 \times 10^{-11}\,\text{m}^3/(\text{kg}\cdot\text{s}^2)$).
• $m_q, m_k$ are the masses of particles $q$ and $k$.
• $\mathbf{s}_q(t), \mathbf{s}_k(t)$ are 2D position vectors.
• $|\mathbf{s}_q(t) - \mathbf{s}_k(t)|$ is the Euclidean distance between particle $k$ and particle $q$:

$|\mathbf{s}_q - \mathbf{s}_k| = \sqrt{(s_{q,X} - s_{k,X})^2 + (s_{q,Y} - s_{k,Y})^2}$

The net gravitational force $\mathbf{F}_q(t)$ acting on particle $q$ is the vector sum of forces from all other $n-1$ particles ($k \in \{0, 1, \dots, n-1\}$ with $k \neq q$):

$$\mathbf{F}_q(t) = \sum_{k \neq q} \mathbf{f}_{qk}(t) = -G m_q \sum_{k \neq q} \dfrac{m_k}{|\mathbf{s}_q(t) - \mathbf{s}_k(t)|^3} \left[ \mathbf{s}_q(t) - \mathbf{s}_k(t) \right]$$

Since acceleration is the second time derivative of position ($\mathbf{a}_q(t) = \mathbf{s}_q''(t)$) and $\mathbf{F}_q = m_q \mathbf{a}_q$, the mass $m_q$ cancels:

$$\mathbf{s}_q''(t) = -G \sum_{j \neq q} \dfrac{m_j}{|\mathbf{s}_q(t) - \mathbf{s}_j(t)|^3} \left[ \mathbf{s}_q(t) - \mathbf{s}_j(t) \right]$$

This represents a coupled system of second-order nonlinear differential equations. Because closed-form analytical solutions exist only for $n = 2$ (the Kepler two-body problem), simulations with $n \ge 3$ require numerical integration.

---

7.1.2 Numerical Integration: Euler’s Method

To advance particle positions and velocities through time, we employ Euler’s method, which approximates a function along its tangent line:

$$g(t_0 + \Delta t) \approx g(t_0) + \Delta t \cdot g'(t_0)$$

flowchart LR
  subgraph Euler["Figure 7.1 and 7.2: Numerical Integration via Euler's Method"]
      direction LR
      T0["Time t = t0<br/>Position s(t0)<br/>Velocity v(t0)"] -->|"Calculate Force F(t0)"| ACC["Acceleration:<br/>a(t0) = F(t0) / m"]
      ACC -->|"Euler Tangent Step"| T1["Time t = t0 + dt<br/>v(t0 + dt) = v(t0) + dt x a(t0)<br/>s(t0 + dt) = s(t0) + dt x v(t0)"]
  end

Applying this linear approximation to velocity and position yields:

$$\mathbf{v}_q(t + \Delta t) \approx \mathbf{v}_q(t) + \Delta t \cdot \mathbf{a}_q(t) = \mathbf{v}_q(t) + \dfrac{\Delta t}{m_q} \mathbf{F}_q(t)$$ $$\mathbf{s}_q(t + \Delta t) \approx \mathbf{s}_q(t) + \Delta t \cdot \mathbf{v}_q(t)$$

At each timestep, the simulation computes all pairwise forces, evaluates accelerations, and updates velocities and positions.

---

7.2 Two Serial Algorithms: Basic vs. Reduced

1. The Basic Algorithm ($O(n^2)$ Force Evaluations)

In the basic approach, each particle $q$ independently iterates through all other particles $k \ne q$ to compute its total force:

for each timestep {
    /* 1. Compute forces */
    for each particle q {
        forces[q][X] = forces[q][Y] = 0.0;
        for each particle k != q {
            x_diff = pos[q][X] - pos[k][X];
            y_diff = pos[q][Y] - pos[k][Y];
            dist = sqrt(x_diff * x_diff + y_diff * y_diff);
            dist_cubed = dist * dist * dist;
            forces[q][X] -= G * masses[q] * masses[k] / dist_cubed * x_diff;
            forces[q][Y] -= G * masses[q] * masses[k] / dist_cubed * y_diff;
        }
    }
    /* 2. Update positions and velocities */
    for each particle q {
        pos[q][X] += delta_t * vel[q][X];
        pos[q][Y] += delta_t * vel[q][Y];
        vel[q][X] += (delta_t / masses[q]) * forces[q][X];
        vel[q][Y] += (delta_t / masses[q]) * forces[q][Y];
    }
}

The force calculation requires $n(n - 1) = n^2 - n$ pairwise force evaluations per timestep.

---

2. The Reduced Algorithm (Halving Arithmetic via Newton’s 3rd Law)

Newton’s third law states that every action generates an equal and opposite reaction:

$\mathbf{f}_{kq} = -\mathbf{f}_{qk}$

This means the $n \times n$ matrix of interparticle forces is skew-symmetric with zeros along the main diagonal:

$$\begin{bmatrix} 0 & \mathbf{f}_{01} & \mathbf{f}_{02} & \cdots & \mathbf{f}_{0,n-1} \\ -\mathbf{f}_{01} & 0 & \mathbf{f}_{12} & \cdots & \mathbf{f}_{1,n-1} \\ -\mathbf{f}_{02} & -\mathbf{f}_{12} & 0 & \cdots & \mathbf{f}_{2,n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -\mathbf{f}_{0,n-1} & -\mathbf{f}_{1,n-1} & -\mathbf{f}_{2,n-1} & \cdots & 0 \end{bmatrix}$$

Rather than evaluating all $n(n-1)$ entries, we need only evaluate the entries strictly above the diagonal ($k > q$), computing exactly:

$\dfrac{n(n-1)}{2} \approx \dfrac{n^2}{2} \text{ force interactions}$

Whenever interaction $\mathbf{f}_{qk}$ is evaluated, it is simultaneously added to $\mathbf{F}_q$ and subtracted from $\mathbf{F}_k$:

/* Program 7.1: Reduced Algorithm for Computing n-Body Forces */
for each particle q
    forces[q] = 0;

for each particle q {
    for each particle k > q {
        x_diff = pos[q][X] - pos[k][X];
        y_diff = pos[q][Y] - pos[k][Y];
        dist = sqrt(x_diff * x_diff + y_diff * y_diff);
        dist_cubed = dist * dist * dist;
        force_qk[X] = G * masses[q] * masses[k] / dist_cubed * x_diff;
        force_qk[Y] = G * masses[q] * masses[k] / dist_cubed * y_diff;

        forces[q][X] += force_qk[X];
        forces[q][Y] += force_qk[Y];
        forces[k][X] -= force_qk[X]; /* Newton's third law */
        forces[k][Y] -= force_qk[Y];
    }
}

Initialization Requirement

Because iteration $q$ updates values in forces[k] for $k > q$, the forces array cannot be zeroed inside the outer loop. It must be zeroed in advance across all particles (e.g., using memset on contiguous vector buffers).

---

7.3 Applying Foster’s Methodology

To parallelize the $n$-body problem, we apply Foster’s four-stage design methodology: Partitioning, Communication, Agglomeration, and Mapping.

flowchart TD
  subgraph Fosters["Foster's Methodology Applied to n-Body Simulation"]
      direction TB
      P["1. Partitioning<br/>Tasks: calculate s(t), v(t), and F(t) for every particle at each timestep"]
      C["2. Communication<br/>Massive data exchange between position updates and force calculations"]
      A["3. Agglomeration<br/>Group all computations for particle q across time into a composite task"]
      M["4. Mapping<br/>Assign particles to cores (Block vs. Cyclic distribution)"]
      P --> C --> A --> M
  end

1. Partitioning and Communication

Initially, primitive tasks compute individual position updates, velocity updates, and force contributions. However, analyzing communication dependencies reveals that calculating $\mathbf{v}_q(t + \Delta t)$ requires $\mathbf{v}_q(t)$ and $\mathbf{F}_q(t)$, and $\mathbf{s}_q(t + \Delta t)$ requires $\mathbf{s}_q(t)$ and $\mathbf{v}_q(t)$.

2. Agglomeration

Because the tightest communication occurs among the variables of a single particle, we agglomerate all calculations associated with particle $q$ ($\mathbf{s}_q, \mathbf{v}_q, \mathbf{F}_q$) into a single composite task.

• In the basic algorithm, particle $q$ receives position $\mathbf{s}_r(t)$ from all other particles $r$.
• In the reduced algorithm, particle $q$ receives position $\mathbf{s}_r(t)$ for $r > q$, and sends the computed interaction $\mathbf{f}_{qr}(t)$ back to particle $r$.

3. Mapping: Block vs. Cyclic Work Partitioning

Since simulation time is sequential, we map spatial particles to cores:

• Basic Algorithm: Every particle performs exactly $n-1$ pairwise checks. A contiguous block partition ($n / \text{thread\_count}$ particles per core) yields perfect load balance and optimal cache spatial locality.
• Reduced Algorithm: Particle $q=0$ performs $n-1$ inner loop iterations, while particle $q=n-1$ performs $0$ iterations! A static block partition causes severe load imbalance (early threads do nearly all the work). A cyclic partition distributes iterations evenly across threads, though at the expense of strided memory accesses.

---

7.4 Parallelizing the Basic Solver with OpenMP

In the basic algorithm, parallelizing the particle loops is straightforward:

#pragma omp parallel num_threads(thread_count)
for (int step = 1; step <= n_steps; step++) {
    #pragma omp single
    if (output_freq) Print_particles(curr, n);

    #pragma omp for schedule(static, n/thread_count)
    for (int q = 0; q < n; q++) {
        Compute_total_force(q, forces, curr, n);
    }

    #pragma omp for schedule(static, n/thread_count)
    for (int q = 0; q < n; q++) {
        Update_particle(q, forces, curr, n, delta_t);
    }
}

Critical Synchronization & Correctness Details:

1. Thread Team Reuse: Placing #pragma omp parallel outside the timestep loop forks the thread team once, avoiding thousands of expensive fork-join cycles.
2. I/O Serialization: The #pragma omp single directive guarantees that exactly one thread executes console output.
3. Loop Barriers: The implicit barrier at the end of #pragma omp for ensures all forces are computed before any thread updates positions, and all positions are updated before the next timestep starts.
4. Data Race Freedom: Thread $q$ only writes to forces[q], pos[q], and vel[q]. Position reads from other particles are read-only during the force calculation phase.

---

7.5 Parallelizing the Reduced Solver with OpenMP

When we attempt to parallelize the reduced algorithm, a major challenge emerges:

$\mathbf{F}_3 = -\mathbf{f}_{03} - \mathbf{f}_{13} - \mathbf{f}_{23}$

If Thread 0 computes $\mathbf{f}_{03}$ and Thread 1 computes $\mathbf{f}_{13}$, both threads attempt to update forces[3] concurrently. This creates a critical data race condition!

Flawed Approaches:

• Single #pragma omp critical: Surrounding force updates with a critical section serializes all force computations across all cores, causing massive contention and running far slower than serial code.
• Array of Particle Locks: Allocating an array of omp_lock_t (one per particle) reduces contention, but lock acquisition/release overhead still degrades performance significantly.

---

The Two-Phase Private Accumulation Strategy

The optimal solution divides force evaluation into two distinct phases:

• Phase 1 (Local Accumulation): Each thread writes its force contributions into its own private subarray loc_forces[my_rank][n][DIM]. Because thread $p$ only writes to row $p$, no synchronization or locking is needed!
• Phase 2 (Global Reduction): An implicit barrier finishes Phase 1. Then, each thread aggregates all partial contributions for its assigned particles.

The implementation in OpenMP is structured as follows:

/* Phase 1: Private accumulation */
#pragma omp for schedule(runtime)
for (int q = 0; q < n; q++) {
    for (int k = q + 1; k < n; k++) {
        Compute_pair_force(q, k, force_qk);
        loc_forces[my_rank][q][X] += force_qk[X];
        loc_forces[my_rank][q][Y] += force_qk[Y];
        loc_forces[my_rank][k][X] -= force_qk[X];
        loc_forces[my_rank][k][Y] -= force_qk[Y];
    }
}

/* Phase 2: Owner aggregation */
#pragma omp for schedule(static, n/thread_count)
for (int q = 0; q < n; q++) {
    forces[q][X] = forces[q][Y] = 0.0;
    for (int t = 0; t < thread_count; t++) {
        forces[q][X] += loc_forces[t][q][X];
        forces[q][Y] += loc_forces[t][q][Y];
    }
}

---

Table 7.1: Phase 1 Computation Grid with Block Partition (3 Threads, 6 Particles)

Thread	Particle	Thread 0 Contributions	Thread 1 Contributions	Thread 2 Contributions
0	0	$f_{01} + f_{02} + f_{03} + f_{04} + f_{05}$	$0$	$0$
0	1	$-f_{01} + f_{12} + f_{13} + f_{14} + f_{15}$	$0$	$0$
1	2	$-f_{02} - f_{12}$	$f_{23} + f_{24} + f_{25}$	$0$
1	3	$-f_{03} - f_{13}$	$-f_{23} + f_{34} + f_{35}$	$0$
2	4	$-f_{04} - f_{14}$	$-f_{24} - f_{34}$	$f_{45}$
2	5	$-f_{05} - f_{15}$	$-f_{25} - f_{35}$	$-f_{45}$

Table 7.2: Phase 1 Computation Grid with Cyclic Partition (3 Threads, 6 Particles)

Thread	Particle	Thread 0 Contributions	Thread 1 Contributions	Thread 2 Contributions
0	0	$f_{01} + f_{02} + f_{03} + f_{04} + f_{05}$	$0$	$0$
1	1	$-f_{01}$	$f_{12} + f_{13} + f_{14} + f_{15}$	$0$
2	2	$-f_{02}$	$-f_{12}$	$f_{23} + f_{24} + f_{25}$
0	3	$-f_{03} + f_{34} + f_{35}$	$-f_{13}$	$-f_{23}$
1	4	$-f_{04} - f_{34}$	$-f_{14} + f_{45}$	$-f_{24}$
2	5	$-f_{05} - f_{35}$	$-f_{15} - f_{45}$	$-f_{25}$

Notice how the cyclic distribution balances the non-zero terms across all three threads much more evenly than the block distribution.

---

OpenMP Benchmark Results

The following table presents runtime results for simulating $n = 400$ particles over $1000$ timesteps on a multicore server:

Table 7.3: OpenMP Runtime Comparison ($n = 400$, $1000$ timesteps)

Threads	Basic Solver	Reduced (Default Sched)	Reduced (Forces Cyclic)	Reduced (All Cyclic)
1	$7.71\,\text{s}$	$3.90\,\text{s}$	$3.90\,\text{s}$	$3.90\,\text{s}$
2	$3.87\,\text{s}$	$2.94\,\text{s}$	$1.98\,\text{s}$	$2.01\,\text{s}$
4	$1.95\,\text{s}$	$1.73\,\text{s}$	$1.01\,\text{s}$	$1.08\,\text{s}$
8	$0.99\,\text{s}$	$0.95\,\text{s}$	$0.54\,\text{s}$	$0.61\,\text{s}$

Key Takeaways:

1. Algorithmic Advantage: The reduced solver is approximately $2\times$ faster than the basic solver across all thread counts due to halving the arithmetic operations.
2. Scheduling Impact: For the reduced solver on 8 threads, using schedule(cyclic) for the force calculation drops execution time from $0.95\,\text{s}$ to $0.54\,\text{s}$ (a $75\%$ improvement) due to superior load balancing.
3. False Sharing Tradeoff: Applying cyclic scheduling to all loops slightly degrades performance ($0.61\,\text{s}$ vs $0.54\,\text{s}$) because strided array writes across threads induce cache coherence invalidations (false sharing).

---

7.6 Parallelizing with POSIX Threads (Pthreads)

Implementing the $n$-body solvers with POSIX Threads follows the same structural logic, with key implementation differences:

• Explicit Loop Scheduling: Without #pragma omp for, threads compute their loop boundaries manually using a helper function:

  void Loop_schedule(int my_rank, int thread_count, int n, int is_cyclic,
                     int* first_p, int* last_p, int* incr_p);

• Explicit Synchronization: Pthreads does not provide implicit loop barriers. Code must invoke pthread_barrier_wait(&barrier) between phases to prevent data races.
• Condition Variable Fallback: If pthread_barrier_t is unavailable on a platform, an equivalent barrier is constructed using a mutex and condition variable:

  void Barrier(void) {
      pthread_mutex_lock(&bar_mut);
      bar_count++;
      if (bar_count == thread_count) {
          bar_count = 0;
          pthread_cond_broadcast(&bar_cond);
      } else {
          while (pthread_cond_wait(&bar_cond, &bar_mut) != 0);
      }
      pthread_mutex_unlock(&bar_mut);
  }