       );
    checkForCudaError("calcA failed");
 
+   threadsInBlock = 256;
+   blockSize = dim3(threadsInBlock);
+   numBlocks = dim3(static_cast<unsigned int>(calcNumBlocks(threadsInBlock, numParticles)));
+   sharedMemoryBytes = d_partLoc.getElementSize() * (threadsInBlock + 1);
+   DevMem<unsigned int, DevMemReuse> d_firstOob(1);
+   d_firstOob.copyArrayToDev(&numParticles, 1);
+   findFirstOob<<<numBlocks, blockSize, sharedMemoryBytes, stream2>>>(
+      sortParticleArray ? d_partLoc.getPtr() : dev_particleLocations.getPtr(),
+      numParticles,
+      OOB_PARTICLE,
+      d_firstOob.getPtr()
+      );
+   checkForCudaError("findFirstOob failed");
+   checkCuda(cudaStreamSynchronize(stream2));
+   unsigned int numRemainingParticles;
+   d_firstOob.copyArrayToHost(&numRemainingParticles, 1);
    if(sortParticleArray)
    {
-      assert(numParticles < 1000000);
-      //std::cout << "Sorting!" << std::endl;
-      threadsInBlock = 256;
-      blockSize = dim3(threadsInBlock);
-      numBlocks = dim3(static_cast<unsigned int>(calcNumBlocks(threadsInBlock, numParticles)));
-      sharedMemoryBytes = d_partLoc.getElementSize() * (threadsInBlock + 1);
-      DevMem<unsigned int, DevMemReuse> d_firstOob(1);
-      d_firstOob.copyArrayToDev(&numParticles, 1);
-      //std::cout << "numParticles=" << numParticles << " OOB_PARTICLE=" << OOB_PARTICLE << std::endl;
-      //std::cout << "blockSize.x=" << blockSize.x << " numBlocks=" << numBlocks.x << " sharedMem=" << sharedMemoryBytes << std::endl;
-      //std::cout << "dev_particleLocations.size()=" << dev_particleLocations.size() << std::endl;
-      findFirstOob<<<numBlocks, blockSize, sharedMemoryBytes, stream2>>>(
-         d_partLoc.getPtr(),
-         numParticles,
-         OOB_PARTICLE,
-         d_firstOob.getPtr()
-         );
-      checkForCudaError("findFirstOob failed");
-      checkCuda(cudaStreamSynchronize(stream2));
-      d_firstOob.copyArrayToHost(&numParticles, 1);
-      assert(numParticles < 1000000);
+      numParticles = numRemainingParticles;
    }
 
 

      dev_a3Sum.getPtr(),
      dev_a4Sum.getPtr(),
      cold,
      numParticles,
      NIJ
      );
   topRowBlockBoundaries.freeMem();

                    static_cast<unsigned int>(calcNumBlocks(threadsY, NY)));
   sharedMemoryBytes = 8 * (threadsX + 1) * sizeof(float);
   uint2 tmpUint2;
   tmpUint2.x = numRemainingParticles;
   tmpUint2.y = 0;
#ifndef NO_THRUST
   DevMem<uint2, ParticleAllocator> topRowBlockBoundaries(numBlocks.x * numBlocks.y, tmpUint2);
