     }
     local_nh.insert(i_ref);
     return local_nh;
+  }
+*/
+
+  __global__ void
+  init_indices_krnl(unsigned int* d_clustering
+                  , unsigned int first_frame_above_threshold) {
+    unsigned int gid = blockIdx.x * blockDim.x + threadIdx.x;
+    if (gid < first_frame_above_threshold) {
+      d_clustering[gid] = gid;
+    }
   }
 
   std::vector<std::size_t>

        }
      }
      // allocate memory and copy sorted coords to GPUs
      for (int i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
        cudaMalloc((void**) &d_coords_sorted[i_gpu]
                 , sizeof(float) * n_rows * n_cols);
        cudaMalloc((void**) &d_local_nh[i_gpu]
                 , sizeof(int) * n_rows);
        cudaMemcpy(d_coords_sorted[i_gpu]
                 , tmp_coords_sorted.data()
                 , sizeof(float) * n_rows * n_cols
                 , cudaMemcpyHostToDevice);
      }
    }
    // indices inside this loop are in order of sorted(!) free energies
    bool neighboring_clusters_merged = false;
    std::set<std::size_t> local_nh;
    while ( ! neighboring_clusters_merged) {
      neighboring_clusters_merged = true;
      logger(std::cout) << "initial merge iteration" << std::endl;
      for (std::size_t i=0; i < first_frame_above_threshold; ++i) {
        if (visited_frames.count(i) == 0) {
          visited_frames.insert(i);
std::cout << i << " / " << first_frame_above_threshold-1 << std::endl;
std::cout << "  high_density_neighborhood" << std::endl;
          // all frames/clusters in local neighborhood should be merged ...
          local_nh = high_density_neighborhood(d_coords_sorted
                                             , n_rows
                                             , n_cols
                                             , d_local_nh
                                             , first_frame_above_threshold
                                             , i
                                             , 4*sigma2
                                             , n_gpus);
          //TODO: profiling!
          //      if this needs lots of time: use OMP parallel sections with
          //      pre-calculated local_nh.
std::cout << "  lump_initial_clusters" << std::endl;
          neighboring_clusters_merged = lump_initial_clusters(local_nh
                                                            , distinct_name
                                                            , clustering
                                                            , fe_sorted
                                                            , first_frame_above_threshold)
                                     && neighboring_clusters_merged;
        }
      }
    }
    // cleanup CUDA environment
    for (int i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaFree(d_coords_sorted[i_gpu]);
      cudaFree(d_local_nh[i_gpu]);
    }
std::cout << "normalize cluster names" << std::endl;
    return normalized_cluster_names(first_frame_above_threshold
                                  , clustering
                                  , fe_sorted);

    // prepare CUDA environment
    int n_gpus = get_num_gpus();
    std::vector<float*> d_coords_sorted(n_gpus);
    std::vector<unsigned int*> d_clustering(n_gpus);
    {
      // sort coords according to free energies
      std::vector<float> tmp_coords_sorted(n_rows * n_cols);
