     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)

               , sizeof(float) * n_rows * n_cols
               , cudaMemcpyHostToDevice);
    }
//    while (microstates_lumped) {



// TODO change state names in clustering to conform to lowest indices

      std::vector<unsigned int> clustering_sorted_orig = clustering_sorted;
      std::cerr << "microstate lumping iteration " << ++i_loop << 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,d_ref_states,\
               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);
        cudaMemcpy(d_ref_states[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]
           , n_rows
           , n_cols
           , max_dist2
           , d_clustering[i_gpu]
           , d_ref_states[i_gpu]
           , i_from
           , i_to);
        }
        cudaDeviceSynchronize();
        check_error("after kernel loop");
      }
      // 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);
      }

      //TODO final lumping

//        for (i=0; i < first_frame_above_threshold; ++i) {
//          if (i_gpu == 0) {
//            clustering_sorted[i] = tmp_clust[i];
//          } else {
//            clustering_sorted[i] = std::min(clustering_sorted[i]
//                                          , tmp_clust[i]);
//          }
//        }


//      std::cout << "lumping CUDA-results ..." << std::endl;
//      // lump microstates
//      clustering_sorted = lumped_clusters(clustering_sorted
//                                        , clustering_sorted_orig
//                                        , first_frame_above_threshold);
//      std::cout << "     ... finished" << std::endl;
//      // compare prev_clustering to clustering_sorted:
//      // if equal, end while
//      microstates_lumped = (clustering_sorted_orig != clustering_sorted);
//    } // end while (if microstates_lumped == false)


    // 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]);
      cudaFree(d_ref_states[i_gpu]);
    }
    // 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];
    }
    return normalized_cluster_names(first_frame_above_threshold
                                  , clustering
                                  , fe_sorted);

    }
    __syncthreads();
    ////////

    // update result for given frame
    if (gid < i_to) {
