 	return d_distance_cache;
 }
 
+static boid * prepare_device_boids(int n) {
+	int boids_bytes = n * sizeof(boid);
+	boid *d_boids = NULL;
+	assert(n > 0);
+	cudaMalloc((void**) &d_boids, boids_bytes);
+	assert(d_boids);
+	return d_boids;
+}
+
 void simulate(simulation_params *sp) {
 	static float *d_distance_cache = NULL;
-	if (!d_distance_cache)
-		d_distance_cache = prepare_distance_cache(sp->boids, sp->n);
-	reload_distance_cache(d_distance_cache, sp->boids, sp->n);
-	calculate_all_forces(sp->boids, sp->n, sp->eps, d_distance_cache);
-	apply_all_forces(sp->boids, sp->n, sp->dt, sp->width, sp->height);
-}
+	static boid *d_boids = NULL;
+	int boids_bytes = sp->n * sizeof(boid);
+	if (!d_distance_cache) {
+		d_distance_cache = prepare_distance_cache(sp->n);
+		d_boids = prepare_device_boids(sp->n);
+	}
+	cudaMemcpy(d_boids, sp->boids, boids_bytes, cudaMemcpyHostToDevice);
+	check_cuda_error();
+	reload_distance_cache(d_distance_cache, d_boids, sp->n);
+	calculate_all_forces(d_boids, sp->n, sp->eps, d_distance_cache);
+	apply_all_forces(d_boids, sp->n, sp->dt, sp->width, sp->height);
+	cudaThreadSynchronize();
+	check_cuda_error();
+	cudaMemcpy(sp->boids, d_boids, boids_bytes, cudaMemcpyDeviceToHost);
+	check_cuda_error();
+}

	}
}

void calculate_all_forces(boid* h_boids, int n, int eps, float *d_distance_cache) {
	int blocksize = 64, neighbours_bytes = n * n * sizeof(int),
		boids_bytes = n * sizeof(boid);
	dim3 threads(blocksize), blocks(n / blocksize);
	static boid *d_boids = NULL;
	static int *d_neighbours = NULL;
	if (!d_boids) {
		cudaMalloc((void**) &d_boids, boids_bytes);
		cudaMalloc((void**) &d_neighbours, neighbours_bytes);
		assert(d_boids && d_neighbours);
	}
	check_cuda_error();
	for (int i = 0; i < n; ++i)
		find_neighbours(d_neighbours + n * i, n, i, d_distance_cache, eps);
	check_cuda_error();
	cudaMemcpy(d_boids, h_boids, boids_bytes, cudaMemcpyHostToDevice);
	check_cuda_error();
	calculate_forces<<<blocks, threads>>>(d_boids,
			d_distance_cache, n, d_neighbours);
	check_cuda_error();
	cudaThreadSynchronize();
	check_cuda_error();
	cudaMemcpy(h_boids, d_boids, boids_bytes, cudaMemcpyDeviceToHost);
	check_cuda_error();
}


		   square(boids[iy].x - boids[ix].x);
}

void reload_distance_cache(float *d_cache, boid *d_boids, int n) {
	const int blocksize = 16;
	dim3 threads(blocksize, blocksize);
	dim3 blocks(n / blocksize, n / blocksize);
	count_distance<<<blocks, threads>>>(d_boids, d_cache, n);
}

__device__ void separation(boid *boids, int self, const int *neighbours,
