 		exit(-1); \
 	} }
 
-__global__ static void neighbourhood(unsigned int *neighbours,
-		unsigned int *flags, const float *d_distances, const int n,
-		const int eps_sq) {
+__global__ static void neighbourhood(int *neighbours, unsigned int *flags,
+		const float *d_distances, const int n, const int eps_sq) {
 	int ix = blockIdx.x * blockDim.x + threadIdx.x,
 		iy = blockIdx.y * blockDim.y + threadIdx.y,
 		offset = ix + iy * n;

	}
}

__global__ static void compact(int *neighbours_out,
		unsigned int *neighbours_in, unsigned int *scanned_flags, int n) {
	const int delim = INT_MAX;
	int ix = blockIdx.x * blockDim.x + threadIdx.x,
		iy = blockIdx.y * blockDim.y + threadIdx.y,
		offset = iy * n;
	if (iy < n) {
		if (0 == ix) {
			if (1 == scanned_flags[offset])
				neighbours_out[offset] = neighbours_in[offset];
			neighbours_out[offset + scanned_flags[offset + n - 1]] = delim;
		} else if (ix < n && scanned_flags[offset + ix] >
				scanned_flags[offset + ix - 1]) {
			neighbours_out[offset + scanned_flags[offset + ix] - 1] =
				neighbours_in[offset + ix];
		}
	}
}

	return planhandle;
}

static void find_neighbours(int *d_neighbours, int n, float *d_distances,
		int eps) {
	static unsigned int *d_flags = NULL;
	const unsigned int blocksize = 16;
	unsigned int flags_bytes = n * n * sizeof(*d_flags);
	static CUDPPHandle planhandle = 0;
	const dim3 threads(blocksize, blocksize);
	dim3 blocks(n / blocksize, n / blocksize);
	if (!d_flags) {
		size_t pitch;
		cudaMallocPitch((void**) &d_flags, &pitch, n * sizeof(*d_flags), n);
		planhandle = prepare_scan_plan(n, pitch / sizeof(*d_flags));
	}
	cudaMemset(d_flags, 0, flags_bytes);
	neighbourhood<<<blocks, threads>>>(d_neighbours, d_flags, d_distances, n,
			eps * eps);
	cudppMultiScan(planhandle, d_flags, d_flags, n, n);
	compact<<<blocks, threads>>>(d_neighbours, d_flags, n);
}

__global__ static void count_distance(float *distance, const int n) {
