GPU n-Body Solvers and Shared-Memory Tiling

On massively parallel GPUs, computing performance is driven by throughput rather than single-thread latency. Modern GPUs feature thousands of arithmetic cores organized across dozens of Streaming Multiprocessors (SMs). To fully saturate the hardware and hide memory access latencies, applications must launch tens or hundreds of thousands of concurrent threads.

In this lesson, we explore how to map the $n$-body problem to NVIDIA GPUs using CUDA, examine why the reduced algorithm’s memory footprint breaks down on GPUs, and introduce Shared-Memory Tiling—a powerful architectural optimization that slashes global DRAM traffic by up to $1024\times$.

---

7.11 Mapping the Basic Solver to CUDA

In MIMD programming (OpenMP and MPI), we assigned large blocks of particles to each core because systems have relatively few processing units ($4$ to $64$). On GPUs, however, hardware schedulers switch between ready warps in single clock cycles to hide latency. Consequently, the optimal mapping strategy assigns one thread per particle:

$\text{Total Threads} = n, \quad \text{Block Count} = \text{blk\_ct} = \dfrac{n}{\text{th\_per\_blk}}$

Each thread computes:

int my_particle = blockIdx.x * blockDim.x + threadIdx.x;
for (int step = 1; step <= n_steps; step++) {
    Compute_force(my_particle, forces, curr, n);
    Update_pos_vel(my_particle, forces, curr, n, delta_t);
}

---

Critical Race Conditions & Grid-Wide Synchronization

The naive single-kernel implementation contains two severe race conditions:

1. Intra-Timestep Race: If Thread A finishes its force calculation and updates $\mathbf{s}_A(t + \Delta t)$ while Thread B is still evaluating forces at time $t$, Thread B will read the future position $\mathbf{s}_A(t + \Delta t)$, corrupting its force calculation.
2. Inter-Timestep Race: If Thread A rushes ahead into timestep $t + 1$ while Thread B is still computing timestep $t$, Thread A reads out-of-date position data for particle B.

Because __syncthreads() only synchronizes threads within the same block (and cannot synchronize across an entire multi-block grid), placing __syncthreads() inside this loop fails.

The Host-Coordinated Grid Barrier Pattern

To enforce a global grid-wide barrier without specialized hardware primitives, we extract the timestep loop into a host function that re-launches distinct kernels:

sequenceDiagram
  autonumber
  participant Host as CPU Host
  participant K1 as Kernel: Compute_force
  participant K2 as Kernel: Update_pos_vel

  loop For each Timestep
      Host->>K1: Launch Compute_force<<<blk_ct, th_per_blk>>>
      Note over K1: All threads compute net forces F[q]
      K1-->>Host: Implicit Grid Barrier (Kernel Completes)
      Host->>K2: Launch Update_pos_vel<<<blk_ct, th_per_blk>>>
      Note over K2: All threads update pos[q] and vel[q]
      K2-->>Host: Implicit Grid Barrier (Kernel Completes)
  end
  Host->>Host: cudaDeviceSynchronize()

/* Host-Coordinated Timestep Driver */
__host__ void Nbody_sim(vect_t forces[], struct particle_s curr[],
                        double delta_t, int n, int n_steps,
                        int blk_ct, int th_per_blk) {
    for (int step = 1; step <= n_steps; step++) {
        /* Kernel 1: Evaluates forces for all particles */
        Compute_force<<<blk_ct, th_per_blk>>>(forces, curr, n);

        /* Implicit barrier: Update_pos_vel waits for Compute_force to finish */
        Update_pos_vel<<<blk_ct, th_per_blk>>>(forces, curr, n, delta_t);
        /* Implicit barrier: Next Compute_force waits for Update_pos_vel */
    }
    cudaDeviceSynchronize();
}

CUDA Cooperative Groups

Starting with CUDA 9 and the NVIDIA Pascal architecture (Compute Capability $\ge 6.0$), developers can synchronize across thread blocks within a single kernel using Cooperative Groups (grid_group::sync()). However, kernel relaunching remains the most portable and universally supported grid barrier technique.

---

7.12 Empirical Performance of the Basic CUDA Solver

The basic CUDA solver was tested on an NVIDIA Pascal GPU (Compute Capability 6.1, 20 SMs, 2560 CUDA cores @ 1.73 GHz) using float precision and $1024$ threads per block:

Table 7.7: Run-times of Basic CUDA vs. Serial CPU n-Body Solver

Blocks (blk_ct)	Particles ($n$)	Serial CPU Run-time	CUDA GPU Run-time	Speedup vs. CPU
1	$1,024$	$0.0701\,\text{s}$	$0.00193\,\text{s}$	$36.3\times$
2	$2,048$	$0.282\,\text{s}$	$0.00289\,\text{s}$	$97.6\times$
32	$32,768$	$50.7\,\text{s}$	$0.0734\,\text{s}$	$690.7\times$
64	$65,536$	$202\,\text{s}$	$0.294\,\text{s}$	$687.1\times$
256	$262,144$	$3,230\,\text{s}$ ($53.8\,\text{min}$)	$2.24\,\text{s}$	$1,442.0\times$
1024	$1,048,576$	$52,000\,\text{s}$ ($14.4\,\text{hours}$)	$48.1\,\text{s}$	$1,081.1\times$

While the serial CPU implementation takes over 14 hours to simulate 1 million particles, CUDA completes the identical workload in just 48 seconds—an astounding $1000\times$ speedup!

---

7.13 Why the Reduced Algorithm Fails on GPUs

On CPU architectures, the reduced algorithm halved execution time by avoiding redundant calculations. Can we port this strategy to CUDA?

Recall that in the reduced algorithm, thread $q$ writes partial force contributions to particle $k$. To avoid race conditions in parallel, each thread requires a private buffer:

$\text{Memory Required} = 2 \times \text{thread\_count} \times n \text{ floats}$

On a GPU where every particle has its own thread ($\text{thread\_count} = n$):

$\text{Total Memory} = 2n^2 \text{ floating-point words}$

For a simulation with $n = 1,000,000$ particles:

$\text{Memory} = 2 \times (10^6)^2 = 2 \times 10^{12} \text{ floats} \times 4\,\text{bytes} = \mathbf{8 \text{ Terabytes of VRAM}!}$

Because high-end GPUs typically provide 16 GB to 80 GB of global memory, the reduced algorithm’s $O(n^2)$ memory footprint is physically impossible to deploy on GPU hardware!

---

7.14 The Solution: Shared-Memory Tiling

To accelerate the basic algorithm without exceeding GPU memory constraints, we optimize memory hierarchy utilization.

In the naive kernel, every thread independently fetches the positions and masses of all $n$ particles from off-chip Global Memory (DRAM). Accessing global DRAM requires 400–800 clock cycles of latency. For $n$ threads, this produces massive memory traffic:

$\text{Global Loads} = 3n(n - 1) \approx 3n^2 \text{ words}$

Instead of relying on hardware L1/L2 caches, we implement Shared-Memory Tiling:

1. Partition the global arrays into logical “tiles” matching the thread block size $\text{tile\_sz} = \text{th\_per\_blk} \le 1024$.
2. In each tile step, the threads of a block collaboratively load $1$ particle each into on-chip __shared__ memory.
3. Synchronize with __syncthreads().
4. All threads evaluate forces against the cached tile data directly from ultra-fast shared memory (1–2 cycles latency).
5. Synchronize with __syncthreads() before loading the next tile.

flowchart LR
  subgraph Tiling["Figure 7.9: Shared-Memory Tile Execution Workflow"]
      direction LR
      GM["Global Memory (DRAM)<br/>Massive Array of n Particles"] -->|"Coalesced Load: 3 words per thread"| SM["Fast On-Chip Shared Memory<br/>Tile t (1024 Particles = 24 KB)"]
      SM -->|"syncthreads Barrier"| COMP["Compute Forces<br/>Inner Loop: 1024 Steps in Fast SRAM"]
      COMP -->|"syncthreads Barrier"| NEXT["Load Next Tile t + 1"]
  end

---

Tiled CUDA Kernel Implementation

#define TH_PER_BLK 1024

__global__ void Compute_force_tiled(
    vect_t forces[],
    const struct particle_s curr[],
    const int n
) {
    __shared__ float s_pos_x[TH_PER_BLK];
    __shared__ float s_pos_y[TH_PER_BLK];
    __shared__ float s_mass[TH_PER_BLK];

    int my_q = blockDim.x * blockIdx.x + threadIdx.x;
    int my_lane = threadIdx.x;

    float my_x = curr[my_q].pos[X];
    float my_y = curr[my_q].pos[Y];
    float my_m = curr[my_q].mass;
    float f_x = 0.0f, f_y = 0.0f;

    int num_tiles = (n + blockDim.x - 1) / blockDim.x;

    for (int t = 0; t < num_tiles; t++) {
        int tile_idx = t * blockDim.x + my_lane;

        /* Collaborative load into on-chip shared memory */
        if (tile_idx < n) {
            s_pos_x[my_lane] = curr[tile_idx].pos[X];
            s_pos_y[my_lane] = curr[tile_idx].pos[Y];
            s_mass[my_lane]  = curr[tile_idx].mass;
        }
        __syncthreads(); /* Wait for tile data to be ready */

        /* Compute interactions against all particles in the tile */
        for (int k = 0; k < blockDim.x; k++) {
            if (t * blockDim.x + k < n && (t * blockDim.x + k) != my_q) {
                float dx = my_x - s_pos_x[k];
                float dy = my_y - s_pos_y[k];
                float dist = sqrtf(dx * dx + dy * dy);
                float dist3 = dist * dist * dist;
                float force = G * my_m * s_mass[k] / dist3;
                f_x -= force * dx;
                f_y -= force * dy;
            }
        }
        __syncthreads(); /* Prevent overwriting tile before all threads finish */
    }

    forces[my_q][X] = f_x;
    forces[my_q][Y] = f_y;
}

---

Memory Bandwidth Reduction Factor

Let’s calculate the reduction in global memory reads:

Basic CUDA Solver: Every thread independently loads $n - 1$ particles from global memory:

$\text{Global Words Loaded}_{\text{Basic}} = 3n(n - 1)$

Tiled CUDA Solver: Each block loads a tile of $\text{th\_per\_blk}$ particles once. Across $\text{blk\_ct}$ tiles, all threads collaboratively load:

$\text{Global Words Loaded}_{\text{Tiled}} = 3n \times \text{blk\_ct} = 3n \left( \dfrac{n}{\text{th\_per\_blk}} \right) = \dfrac{3n^2}{\text{th\_per\_blk}}$

The bandwidth improvement ratio between the two approaches is:

$\text{Bandwidth Improvement Ratio} = \dfrac{3n(n - 1)}{\dfrac{3n^2}{\text{th\_per\_blk}}} \approx \text{th\_per\_blk}$

With $\text{th\_per\_blk} = 1024$, the tiled kernel loads data from global DRAM $1024\times$ fewer times, eliminating memory bus saturation!

---

Empirical Benchmarks: Basic vs. Shared Memory

Both kernels were executed on the same NVIDIA Pascal system with $1024$ threads per block:

Table 7.8: Run-times of Basic vs. Shared-Memory CUDA Solvers

Blocks (blk_ct)	Particles ($n$)	Basic Solver Run-time	Shared-Memory Solver	Speedup Factor
1	$1,024$	$1.93\times 10^{-3}\,\text{s}$	$1.59\times 10^{-3}\,\text{s}$	$1.21\times$
2	$2,048$	$2.89\times 10^{-3}\,\text{s}$	$2.17\times 10^{-3}\,\text{s}$	$1.33\times$
32	$32,768$	$7.34\times 10^{-2}\,\text{s}$	$4.17\times 10^{-2}\,\text{s}$	$1.76\times$
64	$65,536$	$2.94\times 10^{-1}\,\text{s}$	$1.98\times 10^{-1}\,\text{s}$	$1.48\times$
256	$262,144$	$2.24\,\text{s}$	$1.96\,\text{s}$	$1.14\times$
1024	$1,048,576$	$48.1\,\text{s}$	$29.1\,\text{s}$	$1.65\times$

At 1 million particles, the shared-memory tiling optimization reduces execution time from $48.1\,\text{s}$ down to $29.1\,\text{s}$, achieving an additional $1.65\times$ speedup purely by transforming global memory accesses into programmer-managed on-chip cache reuse.