   }
 
   // Save the sum at the appropriate position in the diff vector.
-  if (tix == 0) diff[blockIdx.x*blockDim.y + threadIdx.y] = Z[0][tiy];
+  if (tix == 0) diff[blockIdx.x*blockDim.y + tiy] = Z[0][tiy];
 
 }
 

__global__ void compute_diff(double* X, double* Y, double* diff, int D) {
  
  // Use shared memory for faster computation and reduction
  __shared__ double Z[NUM_THREADS][MAX_THREADS/NUM_THREADS];
  unsigned int tix = threadIdx.x;
  unsigned int tiy = threadIdx.y;
  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;

  if (tix < D) {
    tmp = X[i] - Y[i];
  }

  Z[tix][tiy] = tmp * tmp;


  // Define device arrays and pass data from the CPU to GPU
  double *A, *B, *RetMat;
  CUDA_CHECK(cudaMalloc((void**) &A, D*sizeof(double)));
  CUDA_CHECK(cudaMalloc((void**) &B, Q*D*sizeof(double)));
  CUDA_CHECK(cudaMalloc((void**) &RetMat, Q*NUM_THREADS*sizeof(double)));

  CUDA_CHECK(cudaMemcpy(A, X, D*sizeof(double), 
                        cudaMemcpyHostToDevice));
  CUDA_CHECK(cudaMemcpy(B, Y, Q*D*sizeof(double), 
                        cudaMemcpyHostToDevice));

/*
