title: N-Body CUDA optimization
body: <p>I'm developing an N-body algorithm in CUDA and I would like to learn some tips and tricks for optimization.</p>
  
  <p>I've managed to get <code>16384</code> bodies to run at <code>20Flops</code> on an NVIDIA Geforce GTX 260, which has <code>27</code> Streaming Multiprocessors.</p>
  
  <p>The <code>KernelcomputeForces</code> function is the slow poke taking about <code>95%</code> of the time and I was wondering if there is anymore that I can do to optimize my code.</p>
  
  <p>As far as I see, I've optimized for memory-space-locality and memory-writing. Somewhere in the CUDA docs, it says shared-memory is faster but I dont see how I can make use of that. I've divided the work in <code>16</code> blocks with <code>512</code> threads on each, but thats a bit fuzzy for me.</p>
  
  <p>Please help and thanks for reading this.</p>
  
  <pre><code>n   is number of bodies
  
  gm  is the gpu mass pointer
  
  gpx is the gpu position x pointer
  
  gpy is the gpu position y pointer
  
  gpz is the gpu position z pointer
  
  gfx is the gpu force x pointer
  
  gfy is the gpu force y pointer
  
  gfz is the gpu force z pointer
  </code></pre>
  
  <p><strong>The relevant kernel function</strong></p>
  
  <pre><code>__global__ void KernelcomputeForces( unsigned int n, float* gm, float* gpx, float* gpy, float* gpz, float* gfx, float* gfy, float* gfz ){
      int tid = blockDim.x * blockIdx.x + threadIdx.x;
      int numThreads = blockDim.x * gridDim.x;
  
      float GRAVITY = 0.00001f;
  
      //compare all with all
      for( unsigned int ia=tid; ia&lt;n; ia+=numThreads ){
          float lfx = 0.0f;
          float lfy = 0.0f;
          float lfz = 0.0f;
  
          for( unsigned int ib=0; ib&lt;n; ib++ ){
              //compute distance
              float dx = ( gpx[ib] - gpx[ia]);
              float dy = ( gpy[ib] - gpy[ia] );
              float dz = ( gpz[ib] - gpz[ia] );
              //float distance = sqrt( dx*dx + dy*dy + dz*dz );
              float distanceSquared = dx*dx + dy*dy + dz*dz;
  
              //prevent slingshots and division by zero
              //distance += 0.1f;
              distanceSquared += 0.01f;
  
              //calculate gravitational magnitude between the bodies
              //float magnitude = GRAVITY * ( gm[ia] * gm[ib] ) / ( distance * distance * distance * distance );
              float magnitude = GRAVITY * ( gm[ia] * gm[ib] ) / ( distanceSquared );
  
              //calculate forces for the bodies
              //magnitude times direction
              lfx += magnitude * ( dx );
              lfy += magnitude * ( dy );
              lfz += magnitude * ( dz );
          }
  
          //stores local memory to global memory
          gfx[ia] = lfx;
          gfy[ia] = lfy;
          gfz[ia] = lfz;
      }
  }
  
  extern void GPUcomputeForces( unsigned int n, float* gm, float* gpx, float* gpy, float* gpz, float* gfx, float* gfy, float* gfz ){  
      dim3 gridDim( 16, 1, 1 ); //specifys how many blocks in three possible dimensions
      dim3 blockDim( 512, 1, 1 ); //threads per block
      KernelcomputeForces&lt;&lt;&lt;gridDim, blockDim&gt;&gt;&gt;( n, gm, gpx, gpy, gpz, gfx, gfy, gfz );
  }
  </code></pre>