 		gh_pos[i].z = float((rand())/float(RAND_MAX) - 0.5)*g_Lz;
 		gh_pos[i].w = 0.0f;
 		}
-		
+	
 	// copy particles to the device
 	CUDA_SAFE_CALL(cudaMemcpy(gd_pos, gh_pos, sizeof(float4)*g_N, cudaMemcpyHostToDevice));
 	

	printf("GPU/simple          : ");
	printf("%f ms\n", avg_t);
	}
	
//*************************** simple update method of binning on the GPU
// Run one thread per bin

	// update the data on the device
	cudaMemcpy(gd_pos, gh_pos, sizeof(float4)*g_N, cudaMemcpyHostToDevice);
	}

// sorts the data to mimic HOOMD's standard data pattern (sort of)
void sort_data()
	{
	printf("sorting....\n");
	unsigned int * bin_list = (unsigned int*)malloc(sizeof(unsigned int) *g_N);
	// make even bin dimensions
	float binx = g_Lx / float(g_Mx);
	float biny = g_Ly / float(g_My);
	float binz = g_Lz / float(g_Mz);
	
	float xlo = -g_Lx/2.0f;
	float ylo = -g_Lx/2.0f;
	float zlo = -g_Lx/2.0f;

	// precompute scale factors to eliminate division in inner loop
	float scalex = 1.0f / binx;
	float scaley = 1.0f / biny;
	float scalez = 1.0f / binz;
	
	for (unsigned int i = 0; i < g_N; i++)
		{
		// find the bin each particle belongs in
		unsigned int ib = (unsigned int)((gh_pos[i].x-xlo)*scalex);
		unsigned int jb = (unsigned int)((gh_pos[i].y-ylo)*scaley);
		unsigned int kb = (unsigned int)((gh_pos[i].z-zlo)*scalez);
		
		// need to handle the case where the particle is exactly at the box hi
		if (ib == g_Mx)
			ib = 0;
		if (jb == g_My)
			jb = 0;
		if (kb == g_Mz)
			kb = 0;
			
		// update the bin
		unsigned int bin = ib*(g_Mz*g_My) + jb * g_Mz + kb;
		bin_list[i] = bin;
		}

	bool swapped = false;
	do
		{
		swapped = false;
		for (unsigned int i = 0; i < g_N-1; i++)
			{
			if (bin_list[i] > bin_list[i+1])
				{
				unsigned int tmp = bin_list[i+1];
				bin_list[i+1] = bin_list[i];
				bin_list[i] = tmp;
				
				float4 tmpf = gh_pos[i+1];
				gh_pos[i+1] = gh_pos[i];
				gh_pos[i] = tmpf;
				swapped = true;
				}
			}
		} while (swapped);
	
		
	free(bin_list);
	// update the data on the device
	cudaMemcpy(gd_pos, gh_pos, sizeof(float4)*g_N, cudaMemcpyHostToDevice);
	printf("	done.\n");
	}
	
__global__ void fast_memclear_kernal(unsigned int *d_data, unsigned int N)
	{
	unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
	d_data[idx] = 0;
	}


void rebin_particles_host(unsigned int *idxlist, unsigned int *bin_size, float4 *pos, unsigned int N, float Lx, float Ly, float Lz, unsigned int Mx, unsigned int My, unsigned int Mz, unsigned int Nmax);
	
