 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

HEMI_LAUNCHABLE
void allPairsForcesShared(Vec4f *forceVectors, const Vec4f *bodies, int N)
{
  int idx = blockDim.x * blockIdx.x + threadIdx.x;
  
  extern __shared__ Vec4f jBodies[];

  Vec4f iBody = bodies[idx];
  Vec4f force = Vec4f(0, 0, 0, 0);
 
  for (int tile = 0; tile < gridDim.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

                        const Vec4f *bodies, 
                        int N, bool useShared)
{
  float ss = 0.01f;
  checkCuda( cudaMemcpyToSymbol(HEMI_DEV_CONSTANT(softeningSquared), 
                                &ss, sizeof(float), 0, cudaMemcpyHostToDevice) );

  if (useShared) {
    // we specify the block size and shared memory size, but grid
    // size is automatically chosen
    const int blockSize = 256;
    hemi::ExecutionPolicy ep;
    ep.setBlockSize(blockSize);
    ep.setSharedMemBytes(blockSize * sizeof(Vec4f));
    hemi::cudaLaunch(ep, allPairsForcesShared, forceVectors, bodies, N);
  }
  else // fully automatic configuration in this case
    hemi::cudaLaunch(allPairsForces, forceVectors, bodies, N);
}

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