    subVector(dev_rhoi.getPtr(), dev_rhoe.getPtr(), 
              dev_rho.getPtr(), dev_rhoi.size());
    //#endif
+   first = false;
 }

          int numHotElectrons, int numColdElectrons,
          int numHotIons, int numColdIons)
{

   dim3 *numBlocks;
   dim3 *blockSize;
   int threadsInBlock;

   dev_rhoe.zeroMem();
   dev_rhoi.zeroMem();

   static cudaStream_t stream;
   if(first)
   {
      cudaStreamCreate(&stream);
   }
   // Calculate the rho from the hot electrons
   calcIntermediateRho(dev_rhoe, d_eleHotLoc, 
      numHotElectrons, false, true, stream);
   // Calculate the rho from the cold electrons
   calcIntermediateRho(dev_rhoe, d_eleColdLoc,
      numColdElectrons, true, true, stream);
   // Calculate the rho from the hot ions
   calcIntermediateRho(dev_rhoi, d_ionHotLoc,
      numHotIons, false, false, stream);
   // Calculate the rho from the cold ions
   calcIntermediateRho(dev_rhoi, d_ionColdLoc,
      numColdIons, true, false, stream);

   // Double the rho at the top and bottom of the grid
   threadsInBlock = MAX_THREADS_PER_BLOCK;
