 		V = (int *) malloc(m*n*sizeof(int));
 		out = (int *) malloc(m*k*sizeof(int));
 
-		h_D = (int *) malloc(m*m*sizeof(int));
-//		cudaDeviceSetCacheConfig(cudaFuncCachePreferShared);	
+//h_D = (int *) malloc(m*m*sizeof(int));
+//cudaDeviceSetCacheConfig(cudaFuncCachePreferShared);	
 
 		// compute the execution time
 		cudaEvent_t start, stop;
+//		// create event
+//		cudaEventCreate(&start);
+//		cudaEventCreate(&stop);
+//		// record event
+//		cudaEventRecord(start);
+
+		// allocate space for devices copies
+		cudaMalloc((void **)&d_V, m*n*sizeof(int));
+		cudaMalloc((void **)&d_out, m*k*sizeof(int));
+		cudaMalloc((void **)&D, m*m*sizeof(int));
+
+		for(i=0; i<m*n; i++)
+		{
+			fscanf(fp, "%d", &V[i]);
+		}
+		// copy host values to devices copies
+		cudaMemcpy(d_V, V, m*n*sizeof(int), cudaMemcpyHostToDevice);
+
+		int gridDimX = (int)(ceil((float)m/TILE_WIDTH));
+		int gridDimY = (int)(ceil((float)m/TILE_WIDTH));
+
+		dim3 grid(gridDimX, gridDimY);
+//		dim3 grid(m, m);
+		dim3 block(TILE_WIDTH/2, TILE_WIDTH/2);
+//		dim3 block(TILE_WIDTH, TILE_WIDTH);
+
 		// create event
 		cudaEventCreate(&start);
 		cudaEventCreate(&stop);
 		// record event
 		cudaEventRecord(start);
 
-		// allocate space for devices copies
-		cudaMalloc((void **)&d_V, m*n*sizeof(int));
-		cudaMalloc((void **)&d_out, m*k*sizeof(int));
-		cudaMalloc((void **)&D, m*m*sizeof(int));
-
-		for(i=0; i<m*n; i++)
-		{
-			fscanf(fp, "%d", &V[i]);
-		}
-		// copy host values to devices copies
-		cudaMemcpy(d_V, V, m*n*sizeof(int), cudaMemcpyHostToDevice);
-
-		int gridDimX = (int)(ceil((float)m/TILE_WIDTH));
-		int gridDimY = (int)(ceil((float)m/TILE_WIDTH));
-
-		dim3 grid(gridDimX, gridDimY);
-//		dim3 grid(m, m);
-		dim3 block(TILE_WIDTH/2, TILE_WIDTH/2);
-//		dim3 block(TILE_WIDTH, TILE_WIDTH);
 		// launch knn() kernel on GPU
 		computeDist<<<grid, block>>>(m, n, d_V, D);
 //		computeDist<<<grid, n, n*sizeof(int)>>>(m, n, d_V, D);
 		cudaDeviceSynchronize();
 
-cudaMemcpy(h_D, D, m*m*sizeof(int), cudaMemcpyDeviceToHost);
-//		showD(m, h_D);
+//cudaMemcpy(h_D, D, m*m*sizeof(int), cudaMemcpyDeviceToHost);
+//showD(m, h_D);
 
 		int threadNum = (m<MAX_BLOCK_SIZE)? m: MAX_BLOCK_SIZE;
 //		knn<<<m, m, 2*m*sizeof(int)>>>(m, k, d_V, D, d_out);
-		knn<<<m, threadNum, 2*m*sizeof(int)>>>(m, k, d_V, D, d_out);
-
-		// copy result back to host
-		cudaMemcpy(out, d_out, m*k*sizeof(int), cudaMemcpyDeviceToHost);
+
+		int ptrNumInSMEM = (m<MAX_PTRNUM_IN_SMEM)? m: MAX_PTRNUM_IN_SMEM;
+		knn<<<m, threadNum, 2*ptrNumInSMEM*sizeof(int)>>>(m, k, d_V, D, d_out);
+//		knn<<<m, threadNum, 2*m*sizeof(int)>>>(m, k, d_V, D, d_out);
 
 		// record event and synchronize
 		cudaEventRecord(stop);

/* 
* INPUT:
* m: total num of points
* m is in [10, 1000]
* n: n dimensions
* n is in [1,1000]
* k: num of nearest points
* k is in [1,10]
* V: point coordinates
* the integer elements are in [-5,5]
* OUTPUT:
* out: k nearest neighbors
*/

{
	int i = blockIdx.x;
  	int tid = threadIdx.x;

	int s = blockDim.x/2;
	int resultValue = INIT_MAX;
	int resultIndex = INIT_MAX;
	int indexBase = (m<MAX_PTRNUM_IN_SMEM)? m: MAX_PTRNUM_IN_SMEM;
	
	for(int num=0; num<m; num+=MAX_PTRNUM_IN_SMEM)
	{
		for(int j=tid; j<indexBase; j+=blockDim.x)
		{
			if(j+num == i)
			{
				SMem[j] = INIT_MAX;
			}
			else
			{
				SMem[j] = D[i*m+num+j];
			}
			//index
			SMem[indexBase+j] = j+num;
			__syncthreads();
/*
			if(tid < count)
			{
				if(out[i*k+tid]-num>0 && out[i*k+tid]-num < indexBase)
				{
					SMem[ out[i*k+tid]-num ] = INIT_MAX;
				}
			}
			__syncthreads();
*/
		}
//		__syncthreads();
		for(int j=0; j<count; j++)
		{
			if(out[i*k+j]-num>=0 && out[i*k+j]-num<indexBase)
			{
				SMem[ out[i*k+j]-num ] = INIT_MAX;
			}
			__syncthreads();
		}
		__syncthreads();
/*
		if(tid < count)
		{
			if(out[i*k+tid]-num>0 && out[i*k+tid]-num < indexBase)
			{
				SMem[ out[i*k+tid]-num ] = INIT_MAX;
			}
//			SMem[ out[i*k+tid] ] = INIT_MAX;
		}
		__syncthreads();
*/
		for(s=indexBase/2; s>0; s>>=1) 
		{
			for(int j=tid; j<indexBase; j+=blockDim.x)
			{
				if(j < s) 
				{
					if(SMem[j] == SMem[j+s])
					{
						if(SMem[indexBase+j] > SMem[indexBase+j+s])
						{
							SMem[indexBase+j] = SMem[indexBase+j+s];
						}
					}
					else if(SMem[j] > SMem[j+s])
					{
						SMem[j] = SMem[j+s];
						SMem[indexBase+j] = SMem[indexBase+j+s];
					}
				}
				__syncthreads();
			}
		}
	
		if(resultValue == SMem[0])
		{
			if(resultIndex > SMem[indexBase])
			{
				resultIndex = SMem[indexBase];
			}
		} 
		else if (resultValue > SMem[0])
		{
			resultValue = SMem[0];
			resultIndex = SMem[indexBase];
		}
		__syncthreads();
	}
	return resultIndex;

}

// compute the k nearest neighbors
