Reduce neighbor count computation

cmake/Compadre_Version.txt

@@ -1 +1 @@
-1.7.5
+1.7.6

pycompadre/pycompadre.cpp

@@ -1255,13 +1255,24 @@ Implementation details at:
 
     m_kokkos
     .def("initialize", [](nb::list args) {
+        std::vector<std::string> string_args;
         std::vector<char*> char_args;
+
+        string_args.reserve(args.size());
+        char_args.reserve(args.size() + 1);
+
         for (nb::handle arg : args) {
-            char_args.push_back((char*)nb::cast<std::string>(arg).data());
+            string_args.push_back(nb::cast<std::string>(arg));
+        }
+
+        for (std::string& arg : string_args) {
+            char_args.push_back(const_cast<char*>(arg.c_str()));
         }
+
         char_args.push_back(nullptr);
+
         int argc = (int)args.size();
-        Kokkos::initialize(argc, char_args.data());     
+        Kokkos::initialize(argc, char_args.data());
     })
     .def("finalize", &Kokkos::finalize)
     .def("is_initialized", &Kokkos::is_initialized)

src/Compadre_NeighborLists.hpp

@@ -103,7 +103,6 @@ struct NeighborLists {
 
         Kokkos::deep_copy(_host_cr_neighbor_lists, _cr_neighbor_lists);
         Kokkos::deep_copy(_host_number_of_neighbors_list, _number_of_neighbors_list);
-        Kokkos::deep_copy(_host_row_offsets, _row_offsets);
         Kokkos::fence();
 
         computeRowOffsets();
@@ -133,16 +132,13 @@ struct NeighborLists {
         _host_row_offsets = Kokkos::create_mirror_view(_row_offsets);
 
         Kokkos::deep_copy(_host_number_of_neighbors_list, _number_of_neighbors_list);
-        Kokkos::deep_copy(_host_row_offsets, _row_offsets);
         Kokkos::fence();
 
         computeRowOffsets();
         computeMaxNumNeighbors();
 
         _cr_neighbor_lists = view_type("compressed row neighbor lists data", this->getTotalNeighborsOverAllListsHost());
         _host_cr_neighbor_lists = Kokkos::create_mirror_view(_cr_neighbor_lists);
-        Kokkos::deep_copy(_host_cr_neighbor_lists, _cr_neighbor_lists);
-        Kokkos::fence();
 
         _needs_sync_to_host = false;
     }
@@ -303,6 +299,9 @@ struct NeighborLists {
     //! Get the sum of the number of neighbors of all targets' neighborhoods (device)
     KOKKOS_INLINE_FUNCTION
     global_index_type getTotalNeighborsOverAllListsDevice() const {
+        if (this->getNumberOfTargets()==0) {
+            return 0;
+        }
         return TO_GLOBAL(this->getNumberOfNeighborsDevice(this->getNumberOfTargets()-1)) + this->getRowOffsetDevice(this->getNumberOfTargets()-1);
     }
 ///@}
@@ -323,7 +322,8 @@ NeighborLists<view_type> CreateNeighborLists(view_type neighbor_lists, view_type
 
 //! CreateNeighborLists allows for the construction of an object of type NeighborLists with template deduction
 template <typename view_type>
-NeighborLists<view_type> CreateNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list, view_type neighbor_lists_row_offsets) {
+NeighborLists<view_type> CreateNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list,
+        typename NeighborLists<view_type>::internal_row_offsets_view_type neighbor_lists_row_offsets) {
     return NeighborLists<view_type>(neighbor_lists, number_of_neighbors_list, neighbor_lists_row_offsets);
 }
 

src/Compadre_PointCloudSearch.hpp

@@ -490,7 +490,7 @@ class PointCloudSearch {
                     if (is_dry_run || uniform_radius!=0.0) {
                         number_of_neighbors_list(i) = neighbors_found;
                     } else {
-                        compadre_kernel_assert_debug((neighbors_found==(size_t)number_of_neighbors_list(i)) 
+                        compadre_kernel_assert_release((neighbors_found==(size_t)number_of_neighbors_list(i))
                                 && "Number of neighbors found changed since dry-run.");
                     }
 
@@ -512,8 +512,12 @@ class PointCloudSearch {
             };
             Kokkos::parallel_for("radius search CR", policy, radius_search);
             Kokkos::fence();
-            auto nla = CreateNeighborLists(number_of_neighbors_list);
-            return nla.getTotalNeighborsOverAllListsHost();
+            size_t total_neighbors = 0;
+            for (size_t i=0; i<num_target_sites; ++i) {
+                total_neighbors += static_cast<size_t>(number_of_neighbors_list(i));
+            }
+
+            return total_neighbors;
         }
 
 
@@ -635,7 +639,7 @@ class PointCloudSearch {
                     // the only time the second case using 1e-14 is used is when either zero neighbors or exactly one 
                     // neighbor (neighbor is target site) is found.  when the follow on radius search is conducted, the one
                     // neighbor (target site) will not be found if left at 0, so any positive amount will do, however 1e-14 
-                    // should is small enough to ensure that other neighbors are not found
+                    // should be small enough to ensure that other neighbors are not found
 
                     // needs furthest neighbor's distance for next portion
                     compadre_kernel_assert_release((neighbors_found<neighbor_lists.extent(1) || is_dry_run) 
@@ -784,12 +788,16 @@ class PointCloudSearch {
                     t_min_num_neighbors = (neighbors_found < t_min_num_neighbors) ? neighbors_found : t_min_num_neighbors;
 
                     // scale by epsilon_multiplier to window from location where the last neighbor was found
-                    epsilons(i) = (neighbor_distances(neighbors_found-1) > 0) ?
-                        std::sqrt(neighbor_distances(neighbors_found-1))*epsilon_multiplier : 1e-14*epsilon_multiplier;
-                    // the only time the second case using 1e-14 is used is when either zero neighbors or exactly one 
-                    // neighbor (neighbor is target site) is found.  when the follow on radius search is conducted, the one
-                    // neighbor (target site) will not be found if left at 0, so any positive amount will do, however 1e-14 
-                    // should is small enough to ensure that other neighbors are not found
+                    // Guard the zero-neighbor case before indexing neighbor_distances(neighbors_found-1).
+                    if (neighbors_found > 0 && neighbor_distances(neighbors_found-1) > 0) {
+                        epsilons(i) = std::sqrt(neighbor_distances(neighbors_found-1))*epsilon_multiplier;
+                    } else {
+                        epsilons(i) = 1e-14*epsilon_multiplier;
+                    }
+                    // The 1e-14 case is used when either zero neighbors or exactly one
+                    // neighbor (the target site) is found.  When the follow-on radius search is conducted, the one
+                    // neighbor (target site) will not be found if left at 0, so any positive amount will do, however 1e-14
+                    // should be small enough to ensure that other neighbors are not found
 
                     compadre_kernel_assert_release((epsilons(i)<=max_search_radius || max_search_radius==0 || is_dry_run) 
                             && "max_search_radius given (generally derived from the size of a halo region), \
@@ -813,8 +821,12 @@ class PointCloudSearch {
             generateCRNeighborListsFromRadiusSearch(is_dry_run, trg_pts_view, neighbor_lists, 
                     number_of_neighbors_list, epsilons, 0.0 /*don't set uniform radius*/, max_search_radius);
 
-            auto nla = CreateNeighborLists(number_of_neighbors_list);
-            return nla.getTotalNeighborsOverAllListsHost();
+            size_t total_neighbors = 0;
+            for (size_t i=0; i<num_target_sites; ++i) {
+                total_neighbors += static_cast<size_t>(number_of_neighbors_list(i));
+            }
+
+            return total_neighbors;
         }
 }; // PointCloudSearch