                                          , d_radii2
                                          , n_radii
                                          , d_in_radius[i_stream]);
-      //TODO
-      check_error();
 //TODO debugging
-      cudaDeviceSynchronize();
-      if (i == 0) {
-        std::vector<float> tmp_in_rad(n_radii*n_rows);
-        cudaMemcpy(tmp_in_rad.data()
-                 , d_in_radius[i_stream]
-                 , sizeof(float) * n_rows * n_radii
-                 , cudaMemcpyDeviceToHost);
-        for (auto f: tmp_in_rad) {
-          if (f < 0) {
-            std::cout << "# " << -1.0 * f << std::endl;
-          } else {
-            std::cout << "@  " << f << std::endl;
-          }
-        }
-        exit(EXIT_FAILURE);
-      }
+//      check_error();
+//      cudaDeviceSynchronize();
+//      if (i == 0) {
+//        std::vector<float> tmp_in_rad(n_radii*n_rows);
+//        cudaMemcpy(tmp_in_rad.data()
+//                 , d_in_radius[i_stream]
+//                 , sizeof(float) * n_rows * n_radii
+//                 , cudaMemcpyDeviceToHost);
+//        for (auto f: tmp_in_rad) {
+//          if (f < 0) {
+//            std::cout << "# " << -1.0 * f << std::endl;
+//          } else {
+//            std::cout << "@  " << f << std::endl;
+//          }
+//        }
+//        exit(EXIT_FAILURE);
+//      }
       // compute pops per radius
       for (unsigned int r=0; r < n_radii; ++r) {
         // pops stored col-wise -> just set an offset ...
         offset = r*n_rows;
         //TODO stupid: don't run over all rows for reduction, use boxlimits!
-        reduce_sum<BSIZE> <<< n_rows_ext/BSIZE
+        reduce_sum<BSIZE> <<< Tools::min_multiplicator(n_rows, BSIZE)
                             , BSIZE
                             , 0
                             , streams[i_stream] >>> (offset

      cudaMemset(d_in_radius[i_stream]
               , 0
               , sizeof(float) * n_rows * n_radii);
      // prune range for faster computation
      // (using largest radius in first dimension)
      auto min_max_box = Clustering::Tools::min_max_box(blimits
                                                      , coords[i*n_cols]
                                                      , radii[0]);
      unsigned int offset = min_max_box.first * BSIZE;
      unsigned int rng = (min_max_box.second-min_max_box.first+1);
      in_radius <<< rng
                  , BSIZE
                  , 0



#define BSIZE 128
#define N_STREAMS 2

namespace Clustering {
namespace Density {
