Distributed n-Body Solvers and the Ring Pass Algorithm

Moving from shared-memory systems to distributed-memory clusters requires rethinking how data is stored and exchanged. In distributed memory, no global address space exists: each process has private physical RAM and communicates via explicit message passing across an interconnect network (such as InfiniBand).

In this lesson, we develop both the basic and reduced $n$-body solvers using MPI, explore the elegant Ring Pass algorithm, and analyze why the distributed reduced solver achieves unprecedented memory scalability.

---

7.7 Parallelizing the Basic Solver with MPI

In the basic $n$-body algorithm, each particle $q$ requires the positions and masses of all $n$ particles to compute its net gravitational force $\mathbf{F}_q$. Once forces are known, position and velocity updates are strictly local to particle $q$.

Data Partitioning and Collective Setup

We partition the $n$ particles evenly across $p = \text{comm\_sz}$ processes using a block distribution, where each process owns:

$\text{loc\_n} = \dfrac{n}{p} \text{ particles}$

Because derived datatypes introduce communication overhead compared to basic types, we store particle attributes in separate contiguous 1D arrays rather than an array of structs:

• masses: Stored globally on all processes ($n$ doubles). Since masses never change during Newtonian simulation, Process 0 broadcasts masses once at initialization with MPI_Bcast.
• pos: Process 0 broadcasts the initial positions of all $n$ particles with MPI_Bcast.
• loc_vel: Velocities are strictly local. Process 0 scatters the initial velocity blocks using MPI_Scatter.

---

In-Place Collective Optimization: MPI_IN_PLACE

At each timestep, after processes update their local positions, they must gather updated positions into the global pos array. Normally, an MPI_Allgather call requires separate send and receive buffers:

MPI_Allgather(loc_pos, loc_n, vect_mpi_t, pos, loc_n, vect_mpi_t, comm);

To eliminate the extra buffer allocation for loc_pos, each process can point its local data directly into the corresponding block of the global pos array:

$\text{Process } q \implies \text{loc\_pos} = \text{pos} + q \cdot \text{loc\_n}$

Using the special MPI constant MPI_IN_PLACE, the calling process instructs MPI that its contribution already resides in its receive buffer:

MPI_Allgather(MPI_IN_PLACE, loc_n, vect_mpi_t, pos, loc_n, vect_mpi_t, comm);

---

Optimizing the Outer Timestep Loop

A naive timestep loop might gather positions before printing output, and then gather them again before force calculations. By carefully ordering the operations, we eliminate redundant collective calls:

/* Program 7.2: MPI Basic n-Body Solver Timestep Loop */
Get_input_data(masses, pos, loc_vel, n, loc_n, my_rank, comm);

for (int step = 1; step <= n_steps; step++) {
    if (output_step) {
        Print_particles(pos, loc_vel, n, loc_n, my_rank, comm);
    }

    /* 1. Local force computation using global positions */
    for (int loc_q = 0; loc_q < loc_n; loc_q++) {
        Compute_total_force(loc_q, loc_forces, pos, masses, n, loc_n);
    }

    /* 2. Local Euler integration update */
    for (int loc_q = 0; loc_q < loc_n; loc_q++) {
        Update_particle(loc_q, loc_pos, loc_vel, loc_forces, masses, delta_t);
    }

    /* 3. Global Allgather positioned at the loop end */
    MPI_Allgather(MPI_IN_PLACE, loc_n, vect_mpi_t, pos, loc_n, vect_mpi_t, comm);
}

Placing MPI_Allgather at the end of the timestep loop ensures updated positions are available for output in the next iteration and ready for the subsequent force computation without an extra communication step!

---

7.8 The Reduced Solver and Naive Communication Limits

When parallelizing the reduced solver in MPI, Newton’s third law ($\mathbf{f}_{kq} = -\mathbf{f}_{qk}$) means Process $P_i$ computes forces affecting particles owned by Process $P_j$.

flowchart LR
  subgraph Naive["Figure 7.6: Complex Communication in Naive Reduced MPI"]
      direction LR
      P0["Process 0<br/>Particles 0, 1"] -->|Positions s0, s1| P1["Process 1<br/>Particles 2, 3"]
      P1 -->|Positions s2, s3| P2["Process 2<br/>Particles 4, 5"]
      P0 -->|Forces f02, f03| P1
      P0 -->|Forces f04, f05| P2
      P1 -->|Forces f24, f25| P2
  end

A naive implementation would require:

1. Gathering subsets of positions from higher-ranked processes.
2. Computing partial forces.
3. Scattering partial forces across arbitrary processes and accumulating them into local arrays.

With irregular cyclic distributions, this creates an unmanageable, point-to-point communication graph that degrades throughput and leads to deadlock vulnerabilities.

---

7.9 The Ring Pass Algorithm

The Ring Pass Algorithm solves this communication complexity by mapping processes into a virtual circular ring:

flowchart TD
  subgraph RingTopology["Figure 7.7 and 7.8: Virtual Ring Pipeline"]
      direction TB
      P0["Process 0"] -->|"Send to Lower: (q - 1 + p) mod p"| P3["Process 3"]
      P3 -->|"Send to Lower"| P2["Process 2"]
      P2 -->|"Send to Lower"| P1["Process 1"]
      P1 -->|"Send to Lower"| P0
  end

In a ring topology, each process $q$ only communicates with two immediate neighbors:

• Destination: dest = (q - 1 + comm_sz) % comm_sz (lower-ranked neighbor).
• Source: source = (q + 1) % comm_sz (higher-ranked neighbor).

Algorithmic Concept

Rather than passing only positions, the ring pipeline circulates both positions and forces simultaneously through $p - 1$ communication phases.

In each phase:

1. Each process holds a moving buffer of positions (tmp_pos) and accumulated forces (tmp_forces) received from upstream.
2. The process calculates pairwise interactions between its locally owned particles and the received particles.
3. Computed forces are added to the process’s own loc_forces array and subtracted from the rotating tmp_forces array.
4. The process transmits tmp_pos and tmp_forces downstream to dest and receives the next block from source.
5. After $p$ total phases, a final exchange delivers all remaining forces back to their home processes.

---

Step-by-Step Execution Trace

Consider an example with $n = 4$ particles, $p = 2$ processes, and a cyclic particle distribution:

• Process 0 owns: Particles $0, 2$ ($s_0, s_2$).
• Process 1 owns: Particles $1, 3$ ($s_1, s_3$).

Table 7.4: Ring Pass Force Computation Trace ($n = 4$, $p = 2$)

Phase / Step	Process 0 (loc_forces)	Process 0 (tmp_forces)	Process 1 (loc_forces)	Process 1 (tmp_forces)
Start	$\{0, 0\}$	$\{0, 0\}$	$\{0, 0\}$	$\{0, 0\}$
Local Forces	$\{f_{02}, 0\}$	$\{0, -f_{02}\}$	$\{f_{13}, 0\}$	$\{0, -f_{13}\}$
1st Comm	$\{f_{02}, 0\}$	$\{0, -f_{13}\}$	$\{f_{13}, 0\}$	$\{0, -f_{02}\}$
Compute Phase	$\{f_{01} + f_{02} + f_{03}, f_{23}\}$	$\{-f_{01}, -f_{03} - f_{13} - f_{23}\}$	$\{f_{12} + f_{13}, 0\}$	$\{0, -f_{02} - f_{12}\}$
2nd Comm	$\{f_{01} + f_{02} + f_{03}, f_{23}\}$	$\{0, -f_{02} - f_{12}\}$	$\{f_{12} + f_{13}, 0\}$	$\{-f_{01}, -f_{03} - f_{13} - f_{23}\}$
Final Sum	$\mathbf{F}_0, \mathbf{F}_2$	—	$\mathbf{F}_1, \mathbf{F}_3$	—

Adding tmp_forces into loc_forces at the end completes the net force vector for every particle with zero race conditions!

---

Implementation with MPI_Sendrecv_replace

To minimize message startup latency, tmp_pos and tmp_forces are packed into a single contiguous buffer tmp_data of size $2 \times \text{loc\_n} \times \text{DIM}$ doubles. Communication is executed safely and efficiently using MPI_Sendrecv_replace:

/* Program 7.3: MPI Ring Pass Implementation */
int source = (my_rank + 1) % comm_sz;
int dest   = (my_rank - 1 + comm_sz) % comm_sz;

memcpy(tmp_pos, loc_pos, loc_n * sizeof(vect_t));
memset(loc_forces, 0, loc_n * sizeof(vect_t));
memset(tmp_forces, 0, loc_n * sizeof(vect_t));

/* Step 1: Compute intra-process local interactions */
Compute_local_forces(loc_pos, loc_forces, tmp_forces, loc_n, my_rank, comm_sz);

/* Step 2: Circulate data through the ring */
for (int phase = 1; phase < comm_sz; phase++) {
    MPI_Sendrecv_replace(tmp_data, 2 * loc_n, vect_mpi_t,
                         dest, 0, source, 0, comm, MPI_STATUS_IGNORE);

    int owner = (my_rank + phase) % comm_sz;
    Compute_ring_forces(loc_pos, loc_forces, tmp_pos, tmp_forces,
                        loc_n, my_rank, owner, comm_sz);
}

/* Step 3: Final force hand-off */
MPI_Sendrecv_replace(tmp_forces, loc_n, vect_mpi_t,
                     dest, 0, source, 0, comm, MPI_STATUS_IGNORE);

/* Accumulate received remote contributions */
for (int i = 0; i < loc_n; i++) {
    loc_forces[i][X] += tmp_forces[i][X];
    loc_forces[i][Y] += tmp_forces[i][Y];
}

Global Index Mapping

In a cyclic distribution, local index loc_part1 corresponds to global index glb_part1 = my_rank + loc_part1 * comm_sz. Pairwise forces are evaluated only when glb_part1 < glb_part2, ensuring every interaction is computed exactly once.

---

7.10 Empirical Benchmarks & Memory Scalability Analysis

The MPI solvers were evaluated on an InfiniBand-connected cluster running $n = 800$ particles over $1000$ timesteps:

Table 7.5: MPI n-Body Solver Performance ($n = 800$, $1000$ timesteps)

Cluster Processes	Basic Solver Run-time	Reduced (Ring Pass) Run-time	Reduced Efficiency ($E$)
1	$17.30\,\text{s}$	$8.68\,\text{s}$	$1.00$
2	$8.65\,\text{s}$	$4.45\,\text{s}$	$0.98$
4	$4.35\,\text{s}$	$2.30\,\text{s}$	$0.94$
8	$2.20\,\text{s}$	$1.26\,\text{s}$	$0.86$
16	$1.13\,\text{s}$	$0.78\,\text{s}$	$0.70$

Table 7.6: OpenMP vs. MPI Run-times on 4-Core Node ($n = 800$)

Cores / Threads	OpenMP Basic	OpenMP Reduced	MPI Basic	MPI Reduced
1	$15.13\,\text{s}$	$8.77\,\text{s}$	$17.30\,\text{s}$	$8.68\,\text{s}$
2	$7.62\,\text{s}$	$4.42\,\text{s}$	$8.65\,\text{s}$	$4.45\,\text{s}$
4	$3.85\,\text{s}$	$2.26\,\text{s}$	$4.35\,\text{s}$	$2.30\,\text{s}$

---

The Memory Breakthrough: MPI vs. OpenMP

While OpenMP and MPI achieve comparable runtimes on a single 4-core machine, their memory scaling behaviors diverge drastically as $n$ grows:

• OpenMP Reduced Solver: Requires private force subarrays loc_forces[thread_count][n] to eliminate race conditions, consuming per thread:

$$\text{Memory}_{\text{OpenMP}} = 3\left(\dfrac{n}{p}\right) + 2n \quad \text{[doubles / thread]}$$

• MPI Ring Pass Solver: Each process only stores its local particles plus a single rotating communication buffer of size $4n/p$, consuming per process:

$$\text{Memory}_{\text{MPI}} = n + \dfrac{4n}{p} \quad \text{[doubles / process]}$$

Subtracting the MPI footprint from OpenMP yields the net memory savings per core:

$$\Delta\text{Memory} = \text{Memory}_{\text{OpenMP}} - \text{Memory}_{\text{MPI}} = n - \dfrac{n}{p} \approx n \quad \text{[doubles / core]}$$

For large-scale simulations ($n = 10^7$ particles, $p = 64$ cores):

• OpenMP requires over 10 GB of RAM per thread, quickly exceeding node capacity.
• MPI requires only 160 MB of RAM per process, distributed evenly across cluster nodes.

The MPI Ring Pass algorithm represents a true engineering breakthrough: it halves computational arithmetic while simultaneously enabling simulations orders of magnitude larger than shared-memory systems can physically accommodate.