       cudaFree(d_coords_sorted[i_gpu]);
       cudaFree(d_clustering[i_gpu]);
     }
-std::cout << "normalizing names" << std::endl;
     return normalized_cluster_names(first_frame_above_threshold
                                   , clustering
                                   , fe_sorted);

  }


/*
  __global__ void
  high_density_neighborhood_krnl(std::size_t offset
                               , float* d_coords_sorted
                               , std::size_t first_frame_above_threshold
                               , std::size_t n_cols
                               , std::size_t i_ref
                               , float max_dist2
                               , int* d_local_nh) {
    extern __shared__ float smem[];
    unsigned int bid = blockIdx.x;
    unsigned int tid = threadIdx.x;
    unsigned int bsize = blockDim.x;
    unsigned int gid = bid * bsize + tid + offset;
    // load reference coords into shared memory
    if (tid < n_cols) {
      smem[tid] = d_coords_sorted[i_ref*n_cols+tid];
    }
    // check if in neighborhood
    if (gid < first_frame_above_threshold) {
      float dist2 = 0.0f;
      for (int j=0; j < n_cols; ++j) {
        float c = smem[j] - d_coords_sorted[gid*n_cols+j];
        dist2 = fma(c, c, dist2);
      }
      if (dist2 < max_dist2) {
        d_local_nh[gid] = 1;
      }
    }
  }

  std::set<std::size_t>
  high_density_neighborhood(std::vector<float*> d_coords_sorted
                          , const std::size_t n_rows
                          , const std::size_t n_cols
                          , std::vector<int*> d_local_nh
                          , const std::size_t first_frame_above_threshold
                          , const std::size_t i_ref
                          , const float max_dist2
                          , const int n_gpus) {
    int i_gpu;
    for (i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaMemset(d_local_nh[i_gpu]
               , 0
               , sizeof(int) * n_rows);
    }
    std::size_t rng = first_frame_above_threshold / n_gpus;
    #pragma omp parallel for default(none)\
      private(i_gpu)\
      firstprivate(n_gpus,rng,first_frame_above_threshold,i_ref,max_dist2)\
      shared(d_coords_sorted,d_local_nh)\
      num_threads(n_gpus)\
      schedule(dynamic,1)
    for (i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaSetDevice(i_gpu);
      high_density_neighborhood_krnl
        <<< (rng / BSIZE_SCR) + 1
          , BSIZE_SCR
          , sizeof(float) * n_cols >>>
        (i_gpu * rng
       , d_coords_sorted[i_gpu]
       , first_frame_above_threshold
       , n_cols
       , i_ref
       , max_dist2
       , d_local_nh[i_gpu]);
    }
    cudaDeviceSynchronize();
    // collect results from GPUs
    std::set<std::size_t> local_nh;
    for (i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      std::vector<int> tmp_local_nh(n_rows);
      cudaMemcpy(tmp_local_nh.data()
               , d_local_nh[i_gpu]
               , sizeof(int) * n_rows
               , cudaMemcpyDeviceToHost);
      for (std::size_t i=0; i < n_rows; ++i) {
        if (tmp_local_nh[i] == 1) {
          local_nh.insert(i);
        }
      }
    }
    local_nh.insert(i_ref);
    return local_nh;
  }
*/

  __global__ void
  init_indices_krnl(unsigned int* d_clustering
                  , unsigned int prev_last_frame
                  , unsigned int first_frame_above_threshold) {
    unsigned int gid = blockIdx.x * blockDim.x + threadIdx.x + prev_last_frame;
    if (gid < first_frame_above_threshold) {
      d_clustering[gid] = gid + 1;
    }
  }

  __global__ void
  initial_density_clustering_krnl(unsigned int offset
                                , float* sorted_coords

    __syncthreads();
    if (gid < i_to) {
      unsigned int cluster_id = clustering[gid];
      if (cluster_id == 0) {
        cluster_id = gid + 1;
      }
      unsigned int ref_id = tid+bsize;
      // load reference coordinates for re-use into shared memory
      for (unsigned int j=0; j < n_cols; ++j) {
