 //    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];  
+        } 
+      } 
   }
 
 /*  

                          const int D, const int Q, const int N,
                          double *dist) {

  print_devices_data();

  // Define cuda error
  cudaError_t cudaerr;


  // 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+ MAX_DATA_POINTS*blockIdx.y] 
          = Z[tix][tiy];

}

