 	int *V, *out;				//host copies
 	int *d_V, *d_out;			//device copies
 	int *D;						
+	int *R;						
 //	int *dim_D;					//be replaced with shared memory
 	FILE *fp;
 	if(argc != 2)

		cudaFree(d_V);
		cudaFree(d_out);
		cudaFree(D);
//		cudaFree(dim_D);

		free(V);

			{
				SM[k] += SM[k+s];
			}
			__syncthreads();
		}
	}
	if(k == 0)
	{
		// when n is odd, the last element of SM needs to be added
		if(n > (n/2)*2)
		{
			D[i*m+j] = SM[0] + SM[n-1];
		}
		else
		{
			D[i*m+j] = SM[0];
		}
	}
}

__device__ void initSM(int m, int *D)
{
	int i = blockIdx.x;
	int j = threadIdx.x;
	if(i < j)
	{
		SM[j] = D[i*m+j];
	}
	else
	{
		SM[j] = D[j*m+i];	
	}
}

__device__ int findMin(int m, int n, int k, int count, int *D, int *out, int *R)
{
//	__shared__ int R[1000];
	int i = blockIdx.x;
  	int j = threadIdx.x;
	int s = blockDim.x/2;
	int num;
	int baseIdx = i*m;
	initSM(m, D);
	__syncthreads();
	R[baseIdx+j] = j;
	__syncthreads();
	SM[i] = INIT_MAX;
	__syncthreads();
	for(num=0; num<count; num++)
	{
		SM[ out[i*k+num] ] = INIT_MAX;
	}
	__syncthreads();
	for(s=blockDim.x/2; s>0; s>>=1) 
	{
		// check whether the jth point is the same point as the ith one
		// or has already in the k-nn list
		if(j < s) 
		{
			if(SM[j] > SM[j+s])
			{
				SM[j] = SM[j+s];
				R[baseIdx+j] = R[baseIdx+j+s];
			}
			else if(SM[j] == SM[j+s])
			{
				if(R[baseIdx+j] > R[baseIdx+j+s])
				{
					R[baseIdx+j] = R[baseIdx+j+s];
				}
			}
		}
		__syncthreads();
	}
	// when m is odd, the last element of SM needs to be compared
	if(m > (m/2)*2)
	{
		if(SM[0] > SM[m-1])
		{
			R[baseIdx] = R[baseIdx+m-1];
		}
	}
	__syncthreads();
	return R[baseIdx];
}

// compute the k nearest neighbors
__global__ void knn(int m, int n, int k, int *V, int *D, int *out, int *R)
{
	int i;
	int count;

	i = blockIdx.x;
	for(count=0; count<k; count++)
	{
		out[i*k+count] = findMin(m, n, k, count, D, out, R);
	}
}

