The CUDA Handbook ▸ Chapter 14. N-Body

14.2 Naïve Implementation

Listing 14-1 gives a function that implements the body-body interaction described in the previous section; by annotating it with both the __host__ and __device__ keywords, the CUDA compiler knows it is valid for both the CPU and GPU. The function is templated so it may be invoked for both float and double values (though for this book, only float is fully implemented). It passes back the 3D force vector in the (fx, fy, fz) tuple.

template <typename T>__host__ __device__ void bodyBodyInteraction(    T& ax, T& ay, T& az,    T x0, T y0, T z0,    T x1, T y1, T z1, T mass1,    T softeningSquared){    T dx = x1 - x0;    T dy = y1 - y0;    T dz = z1 - z0;    T distSqr = dx*dx + dy*dy + dz*dz;    distSqr += softeningSquared;    T invDist = (T)1.0 / (T)sqrt(distSqr);    T invDistCube = invDist * invDist * invDist;    T s = mass1 * invDistCube;    ax = dx * s;    ay = dy * s;    az = dz * s;}

Listing 14-2 gives the function that computes the total gravitational force exerted on each body. For each body, it loads that body’s position into (myX, myY, myZ) and then, for every other body, calls bodyBodyInteraction<float> to compute the force exerted between the two. The “AOS” in the function name denotes that the input data comes in the form of an “array of structures”: four packed float values that give the (x, y, z, mass) tuple that specifies a body’s position and mass. The float4 representation is a convenient size for GPU implementation, with native hardware support for loads and stores. Our optimized CPU implementations, described in Section 14.7, make use of so-called “structure of arrays” (SOA) representation where four arrays of float contain packed x, y, z and mass elements for easier processing by SIMD instruction sets. SOA is not a good fit for GPU implementation because the 4 base pointers needed by an SOA representation cost too many registers.

floatComputeGravitation_AOS(    float *force,    float *posMass,    float softeningSquared,    size_t N){    std::chrono::steady_clock::time_point start, end;    start = std::chrono::steady_clock::now();    for ( size_t i = 0; i < N; i++ )    {        float acc[3] = {0, 0, 0};        float myX = posMass[i*4+0];        float myY = posMass[i*4+1];        float myZ = posMass[i*4+2];         for ( size_t j = 0; j < N; j++ ) {            float fx, fy, fz;            float bodyX = posMass[j*4+0];            float bodyY = posMass[j*4+1];            float bodyZ = posMass[j*4+2];            float bodyMass = posMass[j*4+3];             bodyBodyInteraction<float>(                &fx, &fy, &fz,                myX, myY, myZ,                bodyX, bodyY, bodyZ, bodyMass,                softeningSquared );            acc[0] += fx;            acc[1] += fy;            acc[2] += fz;        }         force[3*i+0] = acc[0];        force[3*i+1] = acc[1];        force[3*i+2] = acc[2];    }    end = std::chrono::steady_clock::now();    return (float) std::chrono::duration<double>(end - start).count() * 1000.0f;

Listing 14-3 gives the GPU equivalent to Listing 14-2. For each body, it sums the accelerations due to every other body, then writes that value out to the force array. The L1 and L2 caches on modern GPUs accelerate this workload well, since there is a great deal of reuse in the innermost loop.

Both the outer loop and the inner loop cast the input array posMass to float4, to ensure that the compiler correctly emits a single 16-byte load instruction.

Loop unrolling is an oft-cited optimization for N-Body calculations on GPUs, and it’s not hard to imagine why: branch overhead is much higher on GPUs than CPUs, so the reduced instruction count per loop iteration has a bigger benefit; and the unrolled loop exposes more opportunities for ILP (instruction level parallelism), in which the GPU covers latency of instruction execution as well as memory latency.

To get the benefits of loop unrolling in our N-body application, we need only insert the line:

#pragma unroll <factor>

in front of the for loop over j.

Unfortunately, the optimal loop unrolling factor must be determined empirically. Table 14-1 summarizes the effects of unrolling the loop in this kernel.

Table 14-1. Loop unrolling in the naïve kernel.

In the case of this kernel, in the absence of unrolling, it only delivers 25 billion body-body interactions per second. Even an unroll factor of 2 increases this performance to 30 billion; increasing the unroll factor to 16 delivers the highest performance observed with this kernel: 34.3 billion body-body interactions per second, a 37% performance improvement.

inline __device__ voidComputeNBodyGravitation_Shared_multiGPU(     float *force,     float *posMass,     float softeningSquared,     size_t base,    size_t n,    size_t N ){    extern __shared__ float4 shPosMass[];    for ( int m = blockIdx.x*blockDim.x + threadIdx.x;              m < n;              m += blockDim.x*gridDim.x )    {        size_t i = base+m;        float acc[3] = {0};        float4 myPosMass = ((float4 *) posMass)[i];#pragma unroll 32        for ( int j = 0; j < N; j += blockDim.x ) {            shPosMass[threadIdx.x] = ((float4 *) posMass)[j+threadIdx.x];            __syncthreads();            for ( size_t k = 0; k < blockDim.x; k++ ) {                float fx, fy, fz;                float4 bodyPosMass = shPosMass[k];                 bodyBodyInteraction(                     &fx, &fy, &fz,                     myPosMass.x, myPosMass.y, myPosMass.z,                     bodyPosMass.x,                     bodyPosMass.y,                     bodyPosMass.z,                     bodyPosMass.w,                     softeningSquared );                acc[0] += fx;                acc[1] += fy;                acc[2] += fz;            }            __syncthreads();        }        force[3*m+0] = acc[0];        force[3*m+1] = acc[1];        force[3*m+2] = acc[2];    }}