 #define INIT_MAX 100
 void showResult(int m, int k, int *out);
 
-extern __shared__ int D[];
-
-__device__ void computeDimDist(int m, int n, int *V, int *dim_D)
-{
-	int i = threadIdx.x;
-   	int j = threadIdx.y;
-	int k = threadIdx.z;
-	dim_D[(i*m+j)*n+k] = (V[i*n+k]-V[j*n+k])*(V[i*n+k]-V[j*n+k]);
-}
-
-__device__ void computeDist(int m, int n, int *V, int *dim_D)
-{
-	int i = threadIdx.x;
-   	int j = threadIdx.y;
-	int k = threadIdx.z;
-	int s = blockDim.y/2;
+extern __shared__ int block_dim_D[];
+
+__device__ void computeDimDist(int m, int n, int *V)
+{
+	int i = blockIdx.x;
+   	int j = threadIdx.x;
+	int k = threadIdx.y;
+	block_dim_D[j*n+k] = (V[i*n+k]-V[j*n+k])*(V[i*n+k]-V[j*n+k]);
+}
+
+__device__ void computeDist(int m, int n, int *V, int *D)
+{
+	int i = blockIdx.x;
+   	int j = threadIdx.x;
+	int k = threadIdx.y;
+	int s;
+	int base_idx;
 	int dist = 0;
 	//dimention calculation
-	computeDimDist(m, n, V, dim_D);
+	computeDimDist(m, n, V);
 	__syncthreads();
 	//reduce duplications
-//	if(i < j)
-//	{
-//		for(k=0; k<n; k++)
-//		{
-//			dist += (V[i*n+k]-V[j*n+k])*(V[i*n+k]-V[j*n+k]);
-//		}
-//	}
-	//paralle reduction
+	base_idx = j*n;
 	if(i < j)
 	{
-		for(s=blockDim.y/2; s>0; s>>=1)
+		//paralle reduction
+		for(s=n/2; s>0; s>>=1)
 		{
 			if(k < s)
 			{
-				dim_D[(i*m+j)*n+k] += dim_D[(i*m+j)*n+k+s];
+				block_dim_D[base_idx+k] += block_dim_D[base_idx+k+s];
 			}
 		}
 		__syncthreads();
-		dist += dim_D[(i*m+j)*n];
+		dist = block_dim_D[base_idx];
+		//when n is odd, the last element of dim_D needs to be added
+		if(n > (n/2)*2)
+		{
+			dist += block_dim_D[base_idx+n-1];
+		}
 	}
 	D[i*m+j] = dist;
 }

	int *V, *out;				//host copies
	int *d_V, *d_out;			//device copies
	
//	int *D;						//will be replaced with shared memory
	int *dim_D;

	FILE *fp;
	if(argc != 2)

		cudaFree(d_V);
		cudaFree(d_out);

		cudaFree(D);
//		cudaFree(dim_D);

		free(V);
		free(out);
