N-Body computations are a family of computations that model a set of particles, each of which must consider all the other particles during the computation. Example applications of N-body include (but are not limited to):
gravitational simulation in which stars exert gravitational forces,
molecular modeling in which ions exert electrostatic forces,
particle systems in computer graphics to simulate water and fire, and
“boids,” a technique for computer animation designed to simulate flocking behavior.
Typically the paths of the particles are being simulated per timestep, and computing each timestep costs O(N2) operations for N particles. In most formulations, the forces quickly decrease with distance, leading to hierarchical algorithms in which (for example) the mass and location of the center-of-mass for a collection of particles are used to avoid performing the full O(N2) computations needed otherwise. Barnes-Hut algorithms reduce the runtime to O(NlgN) by introducing a spatial hierarchy that approximates the forces between clusters of objects; for applications where the “leaf nodes” of the computation contain k particles, O(k2) computations must be performed in a given leaf. It is this O(k2) portion of the computation at which GPUs excel.
N-Body workloads have proven the most effective way for GPUs to approach their theoretical limit in processing power. In their GPU Gems 3 paper “Fast N-Body Simulation with CUDA,” Harris et al. frequently cite this theoretical limit in explaining why further performance improvements are not possible. The GPU in question, NVIDIA GeForce 8800 GTX, was so effective at N-body computations that it outperformed custom GRAPE-6 hardware that had been specifically designed to perform astrophysics computation.
In the hopes that readers will be able to “plug in” their computation and find the fastest method for them, this chapter illustrates several different ways to implement N-body and related computations using CUDA:
A naïve implementation illustrates the technique and underscores the effectiveness of caches and the importance of loop unrolling;
A shared memory implementation (for our gravitational computation, the fastest) duplicates Harris et al.’s result, tiling the computation over threadblock-sized collections of bodies to minimize memory latencies in the innermost loop.
A constant memory implementation, inspired by Stone et al.’s implementation of Direct Coulomb Summation (DCS), uses constant memory to hold body descriptions, freeing shared memory for other uses;
A warp-based implementation that exploits the symmetry of forces.
Because readers’ applications may not happen to be gravitational N-body, these different implementations are not presented with the singular goal of optimizing that particular computation. It may make sense to adapt a different implementation depending on the target SM architecture, problem size, and details of the central calculation.
Since gravitational N-body has been presented as a poster child for theoretical performance of GPUs, with speedups of up to 400x reported, the chapter concludes by presenting an implementation optimized for CPUs. By rewriting the calculation to use AVX (Advanced Vector Extensions) and multithreading, a speedup of more than 50x is obtained on an 8-core CPU. Nevertheless, as reported in Section 14.8, the GPU remains roughly an order of magnitude faster than that optimized CPU implementation. In short, the CUDA implementation is faster as well as more readable and maintainable than the optimized CPU implementation.
Most of the chapter’s performance results were measured on a modern desktop machine: an 8-core AMD Ryzen 7 7700X CPU with a GeForce RTX 3060 GPU. The exception is the multi-GPU scaling of Section 14.6, which was measured on the original server-class rig—two Xeon E5-2670 “Sandy Bridge” CPUs with up to four GK104 GPUs, underclocked to conserve power and minimize heat dissipation.
Rather than reporting results in GFLOPS, we report performance results in terms of body-body interactions per second.