 HEMI_LAUNCHABLE
 void allPairsForces(Vec4f *forceVectors, const Vec4f *bodies, int N)
 {
-  int idx = blockDim.x * blockIdx.x + threadIdx.x;
-  forceVectors[idx] = accumulateForce(bodies[idx], bodies, N);
+  for (int idx = blockDim.x * blockIdx.x + threadIdx.x;
+       idx < N;
+       idx += blockDim.x * gridDim.x)
+    forceVectors[idx] = accumulateForce(bodies[idx], bodies, N);
 }
 
 // CUDA kernel that computes all-pairs n-body gravitational forces. Optimized

                        const Vec4f *bodies, 
                        int N, bool useShared)
{
  const int blockSize = 256;
  int gridSize = (N + blockSize - 1) / blockSize;
  hemi::ExecutionPolicy ep;
  ep.setBlockSize(blockSize);
  ep.setGridSize(gridSize);

  float ss = 0.01f;
  checkCuda( cudaMemcpyToSymbol(HEMI_DEV_CONSTANT(softeningSquared), 
                                &ss, sizeof(float), 0, cudaMemcpyHostToDevice) );

  if (useShared) {
    ep.setSharedMemBytes(blockSize * sizeof(Vec4f));
    hemi::cudaLaunch(ep, allPairsForcesShared, forceVectors, bodies, N);
  }
  else
    hemi::cudaLaunch(ep, allPairsForces, forceVectors, bodies, N);
}

// Example of using a host/device class from host code in a .cu file

HEMI_LAUNCHABLE
void allPairsForcesShared(Vec4f *forceVectors, const Vec4f *bodies, int N)
{
  extern __shared__ Vec4f jBodies[];

  for (int idx = blockDim.x * blockIdx.x + threadIdx.x;
       idx < N;
       idx += blockDim.x * gridDim.x)
  {
    Vec4f iBody = bodies[idx];
    Vec4f force = Vec4f(0, 0, 0, 0);
   
    for (int tile = 0; tile < N / blockDim.x; tile++) 
    {
      jBodies[threadIdx.x] = bodies[tile * blockDim.x + threadIdx.x];

      __syncthreads();
          
      force += accumulateForce(iBody, jBodies, blockDim.x);
        
      __syncthreads();
    }

    forceVectors[idx] = force;
  }
}

// Host wrapper function that launches the CUDA kernels
