 __global__ void compute_dist(double* data, double* query, double* dist, 
                               int D, int N) {
   
-  // Use shared memory for faster computation and reduction
+  // Use shared memory for faster computation and reduction (in each block)
   __shared__ double Z[DIM_THREADS][MAX_THREADS/DIM_THREADS];
+  
+  // Thread index in block
   unsigned int tix = threadIdx.x;
   unsigned int tiy = threadIdx.y;
   
-  unsigned int data_ind = threadIdx.x + threadIdx.y*D + blockIdx.x*blockDim.y*D;
+  // Index in current block for each 'data' matrix dimension
+  unsigned int data_ind = threadIdx.x + threadIdx.y*D + 
+                          blockIdx.x*blockDim.y*D;
+
+  // Index in current grid line (y-dimesnion) for every query computation                          
   unsigned int query_ind = threadIdx.x + blockIdx.y*D;
 
-/*  unsigned int i = blockIdx.x*blockDim.x*blockDim.y + threadIdx.y*blockDim.x +
-                    threadIdx.x;
+/*
+  unsigned int i = blockIdx.x*blockDim.x*blockDim.y + 
+                   threadIdx.y*blockDim.x + threadIdx.x;
 */
 
   // Save the difference in the appropriate possition.
   double tmp = 0;
 
+  // When thread's index is less than data's dimenension, then compute distance
   if (tix < D && data_ind < N*D) {
     tmp = data[data_ind] - query[query_ind];
   }
 
   Z[tix][tiy] = tmp * tmp;
 
+  // Synchronize threads before reduction
   __syncthreads();
 
   // Perform reduction

  print_GPU_Mat(deviceQueries, Q*D);
*/  



  compute_dist<<<gridSize, blockSize>>>(deviceData,deviceQueries,
                                          deviceDist,D,N);


                        cudaMemcpyHostToDevice));
  CUDA_CHECK(cudaMemcpy(deviceQueries, queries, Q*D*sizeof(double), 
                        cudaMemcpyHostToDevice));

/*  
  printf("--- Data Matrix ---\n");
