     unsigned int gid = bid * bsize + tid + i_from;
     int comp_size = min(bsize, n_rows - offset);
     if (tid < comp_size) {
-      // reference coordinates to cache
+      // load reference coordinates to cache
       for (unsigned int j=0; j < n_cols; ++j) {
         smem_coords[tid*n_cols+j] = sorted_coords[(tid+offset)*n_cols+j];
       }
-      // reference state information to cache
+      // load reference state information to cache
       smem_cache[col_inter+tid] = clustering[tid+offset];
     }
     __syncthreads();
+
+//TODO: everything lands in state '1'
+
+
     if (gid < i_to) {
-      unsigned int min_id = clustering[gid];
       for (unsigned int j=0; j < n_cols; ++j) {
         // load coordinates of current frame for re-use into shared memory
         smem_coords[(tid+bsize)*n_cols+j] = sorted_coords[gid*n_cols+j];
       }
-      // previous state information to cache
+      // load previous state information to cache
       smem_cache[col_prev+tid] = clustering[gid];
       smem_cache[col_result+tid] = smem_cache[col_prev+tid];
       // compare current frame (tid) against reference block (i)

  }

  std::vector<unsigned int>
  clustering_merged_results(std::vector<std::vector<unsigned int>> clusterings
                          , unsigned int max_row) {
    unsigned int n_results = clusterings.size();
    for (unsigned int i=0; i < max_row; ++i) {
      std::set<unsigned int> start_points;
      for (unsigned int j=0; j < n_results; ++j) {
        start_points.insert(clusterings[j][i]);
      }

      //TODO remove 0 from start_points



      //TODO: finish: follow start_points (id = i_state - 1), collect all
      //              states and rebase them to min(state)
      std::set<unsigned int> update_needed = start_points;
      for (unsigned int s: start_points) {
        unsigned int id = s-1;
        while (clusterings[0][id] != 0
            && clusterings[0][id] != s) {
          update_needed.insert(s);
          s = clusterings[0][id] 

        }


      }
    }
    return clustering;
  }

  std::vector<std::size_t>

               , sizeof(float) * n_rows * n_cols
               , cudaMemcpyHostToDevice);
    }
    // change state names in clustering to conform to lowest indices
    clustering_sorted = clustering_rebased(clustering_sorted);
    // fill zero-set indices with their own index
    for (unsigned int i=0; i < first_frame_above_threshold; ++i) {
      if (clustering_sorted[i] == 0) {
        clustering_sorted[i] = i+1;
      }
    }

//std::cerr << "running CUDA kernels" << std::endl;
    #pragma omp parallel for\
      default(none)\
      private(i,i_gpu,block_rng,i_from,i_to)\
      firstprivate(n_gpus,n_rows,n_cols,gpu_rng,max_dist2,\
                   prev_last_frame,\
                   shared_mem,first_frame_above_threshold)\
      shared(d_coords_sorted,d_clustering,\
             tmp_coords_sorted,clustering_sorted)\
      num_threads(n_gpus)
    for (i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaSetDevice(i_gpu);
      cudaMemcpy(d_clustering[i_gpu]
               , clustering_sorted.data()
               , sizeof(unsigned int) * n_rows
               , cudaMemcpyHostToDevice);
      i_from = prev_last_frame + i_gpu * gpu_rng;
      i_to = (i_gpu == (n_gpus-1))
           ? first_frame_above_threshold
           : prev_last_frame + (i_gpu+1) * gpu_rng;
      block_rng = min_multiplicator(i_to-i_from
                                  , BSIZE_SCR);
      for (i=0; i*BSIZE_SCR < first_frame_above_threshold; ++i) {
        initial_density_clustering_krnl
          <<< block_rng
            , BSIZE_SCR
            , shared_mem >>>
          (i*BSIZE_SCR
         , d_coords_sorted[i_gpu]
         , std::min(first_frame_above_threshold, n_rows)
         , n_cols
         , max_dist2
         , d_clustering[i_gpu]
         , i_from
         , i_to);
      }
      cudaDeviceSynchronize();
      check_error("after kernel loop");
    }
//std::cerr << "collect partial results" << std::endl;
    // collect & merge clustering results from GPUs
    std::vector<std::vector<unsigned int>>
      clstr_results(n_gpus
                  , std::vector<unsigned int>(n_rows));
    for (int i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaMemcpy(clstr_results[i_gpu].data()
               , d_clustering[i_gpu]
               , sizeof(unsigned int) * n_rows
               , cudaMemcpyDeviceToHost);
    }

//for (int i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
//  std::cerr << "partial results from gpu " << i_gpu << std::endl;
//  for (unsigned int i=0; i < first_frame_above_threshold; ++i) {
//    std::cerr << "  " << i << ": " << clstr_results[i_gpu][i] << std::endl;
//  }
//  std::cerr << std::endl;
//}



//std::cerr << "merge partial results" << std::endl;
    clustering_sorted = merge_results(clstr_results
                                    , first_frame_above_threshold);
//std::cerr << "cleanup" << std::endl;
    // cleanup CUDA environment
    for (int i_gpu=0; i_gpu < n_gpus; ++i_gpu) {
      cudaFree(d_coords_sorted[i_gpu]);
      cudaFree(d_clustering[i_gpu]);
    }
//std::cerr << "back-conversion from sorted to true order" << std::endl;
    // convert state trajectory from
    // FE-sorted order to original order
    std::vector<std::size_t> clustering(n_rows);
    for (unsigned int i=0; i < n_rows; ++i) {
      clustering[fe_sorted[i].first] = clustering_sorted[i];
    }
//std::cerr << "normalize and return" << std::endl;
    return normalized_cluster_names(first_frame_above_threshold
                                  , clustering
                                  , fe_sorted);
