 
   // Save the sum at the appropriate position in the diff vector.
   if (tix == 0)// && data_ind < N*D) 
-    dist[threadIdx.y + blockIdx.x*blockDim.y + N*blockIdx.y + offset*MAX_DATA_POINTS] = Z[tix][tiy]; //+ N*blockIdx.y
+    dist[threadIdx.y + blockIdx.x*blockDim.y+ MAX_DATA_POINTS*blockIdx.y] 
+          = Z[tix][tiy];
 
 }
 

  num_blocks_y = Q;
  dim3 gridSize(num_blocks_x,num_blocks_y);


  printf("Block size = (%d,%d), grid size = (%d,%d), D = %d, Q = %d, N = %d\n", 
          DIM_THREADS, block_y_dim, num_blocks_x, num_blocks_y, D, Q, N);
/**/

  // Define and allocate the device space that will hold the appropriate data
  double *deviceData, *deviceQueries, *deviceDist; 
  
  printf("Maximum memory used during computations: %d\n", (D*MAX_DATA_POINTS+Q*D+Q*N)*sizeof(double));

  printf("Allocating device memory for the data matrix.\n");
  CUDA_CHECK(cudaMalloc((void**) &deviceData, 
                        MAX_DATA_POINTS*D*sizeof(double)));
  printf("Allocating device memory for the queries matrix.\n");
  CUDA_CHECK(cudaMalloc((void**) &deviceQueries, Q*D*sizeof(double)));
  printf("Allocating device memory for the distance matrix.\n");
  CUDA_CHECK(cudaMalloc((void**) &deviceDist, Q*N*sizeof(double)));

//  printf("Transfering 'data' matrix from host to device.\n");
//  CUDA_CHECK(cudaMemcpy(deviceData, data, N*D*sizeof(double), 
//                        cudaMemcpyHostToDevice));
  printf("Transfering 'queries' matrix from host to device.\n");
  CUDA_CHECK(cudaMemcpy(deviceQueries, queries, Q*D*sizeof(double), 
                        cudaMemcpyHostToDevice));


//    printf("Call kernel for distance computation.\n");
    compute_dist<<<gridSize, blockSize>>>(deviceData,deviceQueries,
                                          deviceDist,D,N,iter);
    // Check for kernel errors
    cudaerr = cudaGetLastError();
    if (cudaerr != cudaSuccess)
      printf("Error: %s\n", cudaGetErrorString(cudaerr));


//    printf("Transfering 'dist' matrix from device to host.\n");
    CUDA_CHECK(cudaMemcpy(tempDist, deviceDist, Q*MAX_DATA_POINTS*sizeof(double), 
                        cudaMemcpyDeviceToHost));

/*    printf("--- Temp dist ---\n");
    print_gpu_mat(deviceDist, Q*MAX_DATA_POINTS);
    printf("--- Transfer temp dist ---\n");
*/
    for (int qi = 0; qi < Q; qi++) {
        for (int i = 0; i < MAX_DATA_POINTS; i++) {
  //        printf("tempDist[%d] = %f\n", i + qi*MAX_DATA_POINTS, tempDist[i + qi*MAX_DATA_POINTS] );
          dist[i + qi*N + iter*MAX_DATA_POINTS] = tempDist[i + qi*MAX_DATA_POINTS];  
        } 
      } 
  }

/*  
