Stone et al. describe a method of Direct Coulomb Summation (DCS) that
uses shared memory to hold potential map lattice points for a molecular modeling
application, so it must use constant memory to hold body
descriptions. Listing 14-5 shows a CUDA kernel that uses the same method
for our gravitational simulation. Since only 64K of constant memory is
available to developers for a given kernel, each kernel invocation can
only process about 4,000 16-byte body descriptions. The constant
g_bodiesPerPass specifies the number of bodies that can be
considered by the innermost loop.
Since every thread in the innermost loop is reading the same body description, constant memory works well since it is optimized to broadcast reads to all threads in a warp.
const int g_bodiesPerPass = 4000;__constant__ __device__ float4 g_constantBodies[g_bodiesPerPass]; template<typename T>__global__ voidComputeNBodyGravitation_GPU_AOS_const( T *force, T *posMass, T softeningSquared, size_t n, size_t N ){ for ( int i = blockIdx.x*blockDim.x + threadIdx.x; i < N; i += blockDim.x*gridDim.x ) { T acc[3] = {0}; float4 me = ((float4 *) posMass)[i]; T myX = me.x; T myY = me.y; T myZ = me.z; for ( int j = 0; j < n; j++ ) { float4 body = g_constantBodies[j]; float fx, fy, fz; bodyBodyInteraction( &fx, &fy, &fz, myX, myY, myZ, body.x, body.y, body.z, body.w, 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]; }}As shown in Listing 14-6, the host code must loop over the bodies,
calling cudaMemcpyToSymbolAsync() to load the constant memory before
each kernel invocation.
floatComputeNBodyGravitation_GPU_AOS_const( float *force, float *posMass, float softeningSquared, size_t N){ cudaError_t status_cudart; cudaEvent_t evStart = 0, evStop = 0; float ms = 0.0; size_t bodiesLeft = N; void *p; cuda(GetSymbolAddress( &p, g_constantBodies ) ); cuda(EventCreate( &evStart ) ); cuda(EventCreate( &evStop ) ); cuda(EventRecord( evStart, NULL ) ); for ( size_t i = 0; i < N; i += g_bodiesPerPass ) { // bodiesThisPass = max(bodiesLeft, g_bodiesPerPass); size_t bodiesThisPass = bodiesLeft; if ( bodiesThisPass > g_bodiesPerPass ) { bodiesThisPass = g_bodiesPerPass; } cuda(MemcpyToSymbolAsync( g_constantBodies, ((float4 *) posMass)+i, bodiesThisPass*sizeof(float4), 0, cudaMemcpyDeviceToDevice, NULL ) ); ComputeNBodyGravitation_GPU_AOS_const<float> <<<300,256>>>( force, posMass, softeningSquared, bodiesThisPass, N ); bodiesLeft -= bodiesThisPass; } cuda(EventRecord( evStop, NULL ) ); cuda(DeviceSynchronize() ); cuda(EventElapsedTime( &ms, evStart, evStop ) );Error_cudart: cudaEventDestroy( evStop ); cudaEventDestroy( evStart ); return ms;}