              , free_energy.data()
              , sizeof(float) * n_rows
              , cudaMemcpyHostToDevice);
-    // determine optimal block size for nearest neighbor search kernel
     int max_shared_mem;
     cudaDeviceGetAttribute(&max_shared_mem
                          , cudaDevAttrMaxSharedMemoryPerBlock
                          , i_gpu);
     check_error();
-    // max block size B x warp-size (32) given max shared mem M:
-    //   M / (2 * 32 * n_cols * sizeof(float)) = B
-    // -> chosen blocksize = 32*B = M / (2*n_cols*sizeof(float))
-    unsigned int block_size = max_shared_mem / (32*2*n_cols*sizeof(float));
-    // have to do this to get rid of remainder ...
-    block_size *= 32;
-    using Clustering::Tools::min_multiplicator;
+    unsigned int block_size = 128;
+    unsigned int shared_mem = 2 * block_size * n_cols * sizeof(float);
+    if (shared_mem > max_shared_mem) {
+      std::cerr << "error: max. shared mem per block too small on this GPU.\n"
+                << "       either reduce block_size for NN search or get a "
+                <<        "better GPU." << std::endl;
+    }
     unsigned int n_rows_ext = min_multiplicator(n_rows
                                               , block_size)
                             * block_size;
-    unsigned int shared_mem = 2 * block_size * n_cols * sizeof(float);
-    unsigned int rng = min_multiplicator(i_to-i_from, block_size)
-                     * block_size;
-    std::cout << "running kernels"  << std::endl;
-    for (unsigned int i=0; i*block_size < n_rows_ext; ++i) {
+    unsigned int block_rng = min_multiplicator(i_to-i_from, block_size);
+    for (unsigned int i=0; i*block_size <= n_rows_ext; ++i) {
       unsigned int i_stream = i % N_STREAMS_NH;
-
-      std::cout << i_from << " " << i_to << " " << rng << std::endl;
-
-      nearest_neighbor_search <<< rng
+      nearest_neighbor_search <<< block_rng
                                 , block_size
                                 , shared_mem >>> (i*block_size
                                                 , d_coords

                          , std::size_t i_from
                          , std::size_t i_to
                          , int i_gpu) {
    ASSUME_ALIGNED(coords);
    // device buffers
    float* d_coords;

                                                , i_to);
    }
    cudaDeviceSynchronize();
    check_error();
    // initialize neighborhoods
    Neighborhood nh;
    Neighborhood nhhd;
    for (unsigned int i=0; i < n_rows; ++i) {
      nh[i] = {i, std::numeric_limits<float>::max()};
      nhhd[i] = {i, std::numeric_limits<float>::max()};
    }
    // collect partial results from streams
    for (unsigned int i_stream=0; i_stream < N_STREAMS_NH; ++i_stream) {
      std::vector<float> dist_ndx(n_rows * 2);
      auto update_nh = [&dist_ndx,n_rows] (Neighborhood& _nh) -> void {
        for (unsigned int i=0; i < n_rows; ++i) {
          if (dist_ndx[i] < _nh[i].second) {
            _nh[i] = {(unsigned int) dist_ndx[n_rows+i]
                    , dist_ndx[i]};
          }
        }
      };
      cudaMemcpy(dist_ndx.data()
               , d_nh[i_stream]
               , sizeof(float) * n_rows * 2
               , cudaMemcpyDeviceToHost);
      update_nh(nh);
      cudaMemcpy(dist_ndx.data()
               , d_nhhd[i_stream]
               , sizeof(float) * n_rows * 2
               , cudaMemcpyDeviceToHost);
      update_nh(nhhd);
    }
    // device cleanup
    cudaFree(d_coords);
    cudaFree(d_fe);
