Separate calls to algorithms internal kernels (#176)

Author: sbaldu

include/CLUEstering/core/detail/CLUEAlpakaKernels.hpp

@@ -47,22 +47,17 @@ namespace clue {
                                           int32_t point_id) {
       if constexpr (N_ == 0) {
         auto binId = tiles.getGlobalBinByBin(base_vec);
-        // get the size of this bin
         auto binSize = tiles[binId].size();
 
-        // iterate inside this bin
         for (int binIter{}; binIter < binSize; ++binIter) {
           int32_t j{tiles[binId][binIter]};
-          // query N_{dc_}(i)
 
           auto coords_j = getCoords<Ndim>(dev_points, j);
-
           float dist_ij_sq = tiles.distance(coords_i, coords_j);
 
           auto k = kernel(acc, clue::internal::math::sqrt(dist_ij_sq), point_id, j);
           *rho_i += (int)(dist_ij_sq <= dc * dc) * k * dev_points.weight[j];
-
-        }  // end of interate inside this bin
+        }
         return;
       } else {
         for (auto i = search_box[search_box.size() - N_][0];
@@ -87,14 +82,12 @@ namespace clue {
           float rho_i{0.f};
           auto coords_i = getCoords<Ndim>(dev_points, i);
 
-          // Get the extremes of the search box
           clue::SearchBoxExtremes<Ndim> searchbox_extremes;
           for (int dim{}; dim != Ndim; ++dim) {
             searchbox_extremes[dim] =
                 clue::nostd::make_array(coords_i[dim] - dc, coords_i[dim] + dc);
           }
 
-          // Calculate the search box
           clue::SearchBoxBins<Ndim> searchbox_bins;
           dev_tiles.searchBox(searchbox_extremes, searchbox_bins);
 
@@ -121,32 +114,25 @@ namespace clue {
                                                          int32_t point_id) {
       if constexpr (N_ == 0) {
         int binId{tiles.getGlobalBinByBin(base_vec)};
-        // get the size of this bin
         int binSize{tiles[binId].size()};
 
-        // iterate inside this bin
         for (int binIter{}; binIter < binSize; ++binIter) {
           const auto j{tiles[binId][binIter]};
-          // query N'_{dm}(i)
           float rho_j{dev_points.rho[j]};
           bool found_higher{(rho_j > rho_i)};
-          // in the rare case where rho is the same, use detid
           found_higher = found_higher || ((rho_j == rho_i) && (rho_j > 0.f) && (j > point_id));
 
-          // Calculate the distance between the two points
           auto coords_j = getCoords<Ndim>(dev_points, j);
 
           float dist_ij_sq = tiles.distance(coords_i, coords_j);
 
           if (found_higher && dist_ij_sq <= dm_sq) {
-            // find the nearest point within N'_{dm}(i)
             if (dist_ij_sq < *delta_i) {
-              // update delta_i and nearestHigher_i
               *delta_i = dist_ij_sq;
               *nh_i = j;
             }
           }
-        }  // end of interate inside this bin
+        }
 
         return;
       } else {
@@ -183,14 +169,12 @@ namespace clue {
           auto coords_i = getCoords<Ndim>(dev_points, i);
           float rho_i{dev_points.rho[i]};
 
-          // Get the extremes of the search box
           clue::SearchBoxExtremes<Ndim> searchbox_extremes;
           for (int dim{}; dim != Ndim; ++dim) {
             searchbox_extremes[dim] =
                 clue::nostd::make_array(coords_i[dim] - dm, coords_i[dim] + dm);
           }
 
-          // Calculate the search box
           clue::SearchBoxBins<Ndim> searchbox_bins;
           dev_tiles.searchBox(searchbox_extremes, searchbox_bins);
 
@@ -222,13 +206,11 @@ namespace clue {
                                     float rhoc,
                                     int32_t n_points) const {
         for (auto i : alpaka::uniformElements(acc, n_points)) {
-          // initialize cluster_index
           dev_points.cluster_index[i] = -1;
 
           float delta_i = dev_points.delta[i];
           float rho_i = dev_points.rho[i];
 
-          // Determine whether the point is a seed or an outlier
           bool is_seed = (delta_i > seed_dc) && (rho_i >= rhoc);
 
           if (is_seed) {
@@ -255,25 +237,18 @@ namespace clue {
 
           int idx_this_seed = seeds_0[idx_cls];
           dev_points.cluster_index[idx_this_seed] = idx_cls;
-          // push_back idThisSeed to localStack
           local_stack[local_stack_size] = idx_this_seed;
           ++local_stack_size;
-          // process all elements in localStack
           while (local_stack_size > 0) {
-            // get last element of localStack
             int idx_end_of_local_stack{local_stack[local_stack_size - 1]};
             int temp_cluster_index = dev_points.cluster_index[idx_end_of_local_stack];
-            // pop_back last element of localStack
             local_stack[local_stack_size - 1] = -1;
             --local_stack_size;
             const auto& followers_ies = followers[idx_end_of_local_stack];
             const auto followers_size = followers_ies.size();
-            // loop over followers of last element of localStack
             for (int j{}; j != followers_size; ++j) {
-              // pass id to follower
               int follower = followers_ies[j];
               dev_points.cluster_index[follower] = temp_cluster_index;
-              // push_back follower to localStack
               local_stack[local_stack_size] = follower;
               ++local_stack_size;
             }
@@ -282,5 +257,60 @@ namespace clue {
       }
     };
 
+    using WorkDiv = clue::WorkDiv<clue::Dim1D>;
+
+    template <concepts::accelerator TAcc, concepts::queue TQueue, uint8_t Ndim, typename KernelType>
+    inline void computeLocalDensity(TQueue& queue,
+                                    const WorkDiv& work_division,
+                                    TilesAlpakaView<Ndim>& tiles,
+                                    PointsView& dev_points,
+                                    KernelType&& kernel,
+                                    float dc,
+                                    int32_t size) {
+      alpaka::exec<TAcc>(queue,
+                         work_division,
+                         KernelCalculateLocalDensity{},
+                         tiles,
+                         dev_points,
+                         std::forward<KernelType>(kernel),
+                         dc,
+                         size);
+    }
+
+    template <concepts::accelerator TAcc, concepts::queue TQueue, uint8_t Ndim>
+    inline void computeNearestHighers(TQueue& queue,
+                                      const WorkDiv& work_division,
+                                      TilesAlpakaView<Ndim>& tiles,
+                                      PointsView& dev_points,
+                                      float dm,
+                                      int32_t size) {
+      alpaka::exec<TAcc>(
+          queue, work_division, KernelCalculateNearestHigher{}, tiles, dev_points, dm, size);
+    }
+
+    template <concepts::accelerator TAcc, concepts::queue TQueue>
+    inline void findClusterSeeds(TQueue& queue,
+                                 const WorkDiv& work_division,
+                                 VecArray<int32_t, reserve>* seeds,
+                                 PointsView& dev_points,
+                                 float seed_dc,
+                                 float rhoc,
+                                 int32_t size) {
+      alpaka::exec<TAcc>(
+          queue, work_division, KernelFindClusters{}, seeds, dev_points, seed_dc, rhoc, size);
+    }
+
+    template <concepts::accelerator TAcc, concepts::queue TQueue>
+    inline void assignPointsToClusters(TQueue& queue,
+                                       std::size_t block_size,
+                                       VecArray<int32_t, reserve>* seeds,
+                                       clue::FollowersView followers,
+                                       PointsView dev_points) {
+      const Idx grid_size = clue::divide_up_by(reserve, block_size);
+      const auto work_division = clue::make_workdiv<TAcc>(grid_size, block_size);
+      alpaka::exec<TAcc>(
+          queue, work_division, KernelAssignClusters{}, seeds, followers, dev_points);
+    }
+
   }  // namespace detail
 }  // namespace clue

include/CLUEstering/core/detail/Clusterer.hpp

@@ -315,51 +315,25 @@ namespace clue {
                                            const KernelType& kernel,
                                            Queue& queue,
                                            std::size_t block_size) {
-    const auto nPoints = h_points.size();
-    // fill the tiles
-    m_tiles->template fill<Acc>(queue, dev_points, nPoints);
-
-    const Idx grid_size = clue::divide_up_by(nPoints, block_size);
-    auto working_div = clue::make_workdiv<Acc>(grid_size, block_size);
-    alpaka::exec<Acc>(queue,
-                      working_div,
-                      detail::KernelCalculateLocalDensity{},
-                      m_tiles->view(),
-                      dev_points.view(),
-                      kernel,
-                      m_dc,
-                      nPoints);
-    alpaka::exec<Acc>(queue,
-                      working_div,
-                      detail::KernelCalculateNearestHigher{},
-                      m_tiles->view(),
-                      dev_points.view(),
-                      m_dm,
-                      nPoints);
+    const auto n_points = h_points.size();
+    m_tiles->template fill<Acc>(queue, dev_points, n_points);
+
+    const Idx grid_size = clue::divide_up_by(n_points, block_size);
+    auto work_division = clue::make_workdiv<Acc>(grid_size, block_size);
+
+    detail::computeLocalDensity<Acc>(
+        queue, work_division, m_tiles->view(), dev_points.view(), kernel, m_dc, n_points);
+    detail::computeNearestHighers<Acc>(
+        queue, work_division, m_tiles->view(), dev_points.view(), m_dm, n_points);
 
     m_followers->template fill<Acc>(queue, dev_points);
 
-    alpaka::exec<Acc>(queue,
-                      working_div,
-                      detail::KernelFindClusters{},
-                      m_seeds->data(),
-                      dev_points.view(),
-                      m_seed_dc,
-                      m_rhoc,
-                      nPoints);
-
-    // We change the working division when assigning the clusters
-    const Idx grid_size_seeds = clue::divide_up_by(reserve, block_size);
-    auto working_div_seeds = clue::make_workdiv<Acc>(grid_size_seeds, block_size);
-
-    alpaka::exec<Acc>(queue,
-                      working_div_seeds,
-                      detail::KernelAssignClusters{},
-                      m_seeds->data(),
-                      m_followers->view(),
-                      dev_points.view());
-    alpaka::wait(queue);
+    detail::findClusterSeeds<Acc>(
+        queue, work_division, m_seeds->data(), dev_points.view(), m_seed_dc, m_rhoc, n_points);
+    detail::assignPointsToClusters<Acc>(
+        queue, block_size, m_seeds->data(), m_followers->view(), dev_points.view());
 
+    alpaka::wait(queue);
     clue::copyToHost(queue, h_points, dev_points);
   }
 
@@ -369,48 +343,23 @@ namespace clue {
                                            const KernelType& kernel,
                                            Queue& queue,
                                            std::size_t block_size) {
-    const auto nPoints = dev_points.size();
-    m_tiles->template fill<Acc>(queue, dev_points, nPoints);
-
-    const Idx grid_size = clue::divide_up_by(nPoints, block_size);
-    auto working_div = clue::make_workdiv<Acc>(grid_size, block_size);
-    alpaka::exec<Acc>(queue,
-                      working_div,
-                      detail::KernelCalculateLocalDensity{},
-                      m_tiles->view(),
-                      dev_points.view(),
-                      kernel,
-                      m_dc,
-                      nPoints);
-    alpaka::exec<Acc>(queue,
-                      working_div,
-                      detail::KernelCalculateNearestHigher{},
-                      m_tiles->view(),
-                      dev_points.view(),
-                      m_dm,
-                      nPoints);
+    const auto n_points = dev_points.size();
+    m_tiles->template fill<Acc>(queue, dev_points, n_points);
+
+    const Idx grid_size = clue::divide_up_by(n_points, block_size);
+    auto work_division = clue::make_workdiv<Acc>(grid_size, block_size);
+
+    detail::computeLocalDensity<Acc>(
+        queue, work_division, m_tiles->view(), dev_points.view(), kernel, m_dc, n_points);
+    detail::computeNearestHighers<Acc>(
+        queue, work_division, m_tiles->view(), dev_points.view(), m_dm, n_points);
 
     m_followers->template fill<Acc>(queue, dev_points);
 
-    alpaka::exec<Acc>(queue,
-                      working_div,
-                      detail::KernelFindClusters{},
-                      m_seeds->data(),
-                      dev_points.view(),
-                      m_seed_dc,
-                      m_rhoc,
-                      nPoints);
-
-    // We change the working division when assigning the clusters
-    const Idx grid_size_seeds = clue::divide_up_by(reserve, block_size);
-    auto working_div_seeds = clue::make_workdiv<Acc>(grid_size_seeds, block_size);
-
-    alpaka::exec<Acc>(queue,
-                      working_div_seeds,
-                      detail::KernelAssignClusters{},
-                      m_seeds->data(),
-                      m_followers->view(),
-                      dev_points.view());
+    detail::findClusterSeeds<Acc>(
+        queue, work_division, m_seeds->data(), dev_points.view(), m_seed_dc, m_rhoc, n_points);
+    detail::assignPointsToClusters<Acc>(
+        queue, block_size, m_seeds->data(), m_followers->view(), dev_points.view());
   }
 
 }  // namespace clue