dynamic backend implementation

CMakeLists.txt

@@ -28,10 +28,6 @@ set(GLM_INC_DIR ${PROJECT_SOURCE_DIR}/third_party/glm)
 
 set_property(CACHE THRUST_BACKEND PROPERTY STRINGS CUDA OMP CPP)
 
-if(NOT THRUST_BACKEND STREQUAL "CUDA" AND NOT THRUST_BACKEND STREQUAL "OMP" AND NOT THRUST_BACKEND STREQUAL "CPP")
-    message(FATAL_ERROR "Invalid value for THRUST_BACKEND: ${THRUST_BACKEND}. Should be one of \"CUDA\", \"OMP\" or \"CPP\"")
-endif()
-
 if(MSVC)
     set(MANIFOLD_FLAGS)
 else()
@@ -41,21 +37,30 @@ set(MANIFOLD_NVCC_FLAGS -Xcudafe --diag_suppress=esa_on_defaulted_function_ignor
 set(MANIFOLD_NVCC_RELEASE_FLAGS -O3)
 set(MANIFOLD_NVCC_DEBUG_FLAGS -G)
 
-if(THRUST_BACKEND STREQUAL "CUDA")
+set(THRUST_INC_DIR ${PROJECT_SOURCE_DIR}/third_party/thrust)
+
+message("Using ${THRUST_BACKEND} backend")
+if (THRUST_BACKEND STREQUAL "CUDA")
     # we cannot set THRUST_INC_DIR when building with CUDA, otherwise the
     # compiler will not use the system CUDA headers which causes incompatibility
     enable_language(CUDA)
-else()
-    set(THRUST_INC_DIR ${PROJECT_SOURCE_DIR}/third_party/thrust)
-    if(THRUST_BACKEND STREQUAL "OMP")
-        message("------------------------- Using OpenMP instead of CUDA.")
-        find_package(OpenMP REQUIRED)
-        set(MANIFOLD_OMP_INCLUDE OpenMP::OpenMP_CXX)
-        set(MANIFOLD_FLAGS ${MANIFOLD_FLAGS} -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_OMP -fopenmp)
-    else()
-        message("------------------------- Using C++ instead of CUDA.")
-        set(MANIFOLD_FLAGS ${MANIFOLD_FLAGS} -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_CPP)
-    endif()
+    find_package(CUDA  REQUIRED)
+    # clear THRUST_INC_DIR, we use the one from nvcc
+    set(THRUST_INC_DIR "")
+    # OpenMP also required
+    find_package(OpenMP REQUIRED)
+    set(MANIFOLD_INCLUDE OpenMP::OpenMP_CXX)
+    set(MANIFOLD_FLAGS ${MANIFOLD_FLAGS} -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_CUDA)
+    set(MANIFOLD_NVCC_FLAGS ${MANIFOLD_NVCC_FLAGS} -fopenmp)
+elseif (THRUST_BACKEND STREQUAL "OMP")
+    find_package(OpenMP REQUIRED)
+    set(MANIFOLD_INCLUDE OpenMP::OpenMP_CXX)
+    set(MANIFOLD_FLAGS ${MANIFOLD_FLAGS} -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_OMP -fopenmp)
+elseif (THRUST_BACKEND STREQUAL "CPP")
+    set(MANIFOLD_FLAGS ${MANIFOLD_FLAGS} -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_CPP)
+else ()
+    message(FATAL_ERROR "Invalid value for THRUST_BACKEND: ${THRUST_BACKEND}. "
+    "Should be one of \"CUDA\", \"OMP\", \"CPP\"")
 endif()
 
 add_subdirectory(utilities)

collider/CMakeLists.txt

@@ -22,7 +22,7 @@ target_include_directories( ${PROJECT_NAME}
 )
 target_link_libraries( ${PROJECT_NAME}
     PUBLIC utilities
-    PRIVATE ${MANIFOLD_OMP_INCLUDE}
+    PRIVATE ${MANIFOLD_INCLUDE}
 )
 
 target_compile_features(${PROJECT_NAME} 

collider/src/collider.cpp

@@ -14,6 +14,7 @@
 
 #include "collider.h"
 #include "utils.h"
+#include "par.h"
 
 #ifdef _MSC_VER
 #include <intrin.h>
@@ -252,7 +253,7 @@ Collider::Collider(const VecDH<Box>& leafBB,
   nodeParent_.resize(num_nodes, -1);
   internalChildren_.resize(leafBB.size() - 1, thrust::make_pair(-1, -1));
   // organize tree
-  thrust::for_each_n(thrust::device, countAt(0), NumInternal(),
+  for_each_n(autoPolicy(NumInternal()), countAt(0), NumInternal(),
                      CreateRadixTree({nodeParent_.ptrD(),
                                       internalChildren_.ptrD(), leafMorton}));
   UpdateBoxes(leafBB);
@@ -273,8 +274,8 @@ SparseIndices Collider::Collisions(const VecDH<T>& querriesIn) const {
     // scalar number of overlaps found
     VecDH<int> nOverlapsD(1, 0);
     // calculate Bounding Box overlaps
-    thrust::for_each_n(
-        thrust::device, zip(querriesIn.cbegin(), countAt(0)), querriesIn.size(),
+    for_each_n(
+        autoPolicy(querriesIn.size()), zip(querriesIn.cbegin(), countAt(0)), querriesIn.size(),
         FindCollisions<T>({querryTri.ptrDpq(), nOverlapsD.ptrD(), maxOverlaps,
                            nodeBBox_.ptrD(), internalChildren_.ptrD()}));
     nOverlaps = nOverlapsD[0];
@@ -305,13 +306,13 @@ void Collider::UpdateBoxes(const VecDH<Box>& leafBB) {
   // copy in leaf node Boxs
   strided_range<VecDH<Box>::Iter> leaves(nodeBBox_.begin(), nodeBBox_.end(),
                                           2);
-  thrust::copy(thrust::device, leafBB.cbegin(), leafBB.cend(), leaves.begin());
+  auto policy = autoPolicy(NumInternal());
+  copy(policy, leafBB.cbegin(), leafBB.cend(), leaves.begin());
   // create global counters
-  VecDH<int> counter(NumInternal());
-  thrust::fill(thrust::device, counter.begin(), counter.end(), 0);
+  VecDH<int> counter(NumInternal(), 0);
   // kernel over leaves to save internal Boxs
-  thrust::for_each_n(
-      thrust::device, countAt(0), NumLeaves(),
+  for_each_n(
+      policy, countAt(0), NumLeaves(),
       BuildInternalBoxes({nodeBBox_.ptrD(), counter.ptrD(), nodeParent_.ptrD(),
                           internalChildren_.ptrD()}));
 }
@@ -330,7 +331,7 @@ bool Collider::Transform(glm::mat4x3 transform) {
     if (count != 2) axisAligned = false;
   }
   if (axisAligned) {
-    thrust::for_each(thrust::device, nodeBBox_.begin(), nodeBBox_.end(),
+    for_each(autoPolicy(nodeBBox_.size()), nodeBBox_.begin(), nodeBBox_.end(),
                      TransformBox({transform}));
   }
   return axisAligned;

flake.nix

@@ -48,10 +48,10 @@
           packages.manifold-omp = manifold { backend = "OMP"; runtime = [ pkgs.llvmPackages.openmp ]; };
           packages.manifold-cuda = manifold {
             backend = "CUDA";
-            runtime = [
-              pkgs.cudaPackages.cudatoolkit_11
+            runtime = with pkgs; [
+              llvmPackages.openmp
+              cudaPackages.cudatoolkit_11
             ];
-            doCheck = false;
           };
           packages.manifold-js = pkgs.buildEmscriptenPackage {
             name = "manifold-js";

manifold/CMakeLists.txt

@@ -30,7 +30,7 @@ endif()
 target_include_directories( ${PROJECT_NAME} PUBLIC ${PROJECT_SOURCE_DIR}/include )
 target_link_libraries( ${PROJECT_NAME}
     PUBLIC utilities
-    PRIVATE collider polygon ${MANIFOLD_OMP_INCLUDE} graphlite
+    PRIVATE collider polygon ${MANIFOLD_INCLUDE} graphlite
 )
 
 target_compile_features(${PROJECT_NAME} 

manifold/src/boolean3.cpp

@@ -13,6 +13,7 @@
 // limitations under the License.
 
 #include "boolean3.h"
+#include "par.h"
 #include <limits>
 
 // TODO: make this runtime configurable for quicker debug
@@ -85,12 +86,13 @@ struct CopyFaceEdges {
 SparseIndices Filter11(const Manifold::Impl &inP, const Manifold::Impl &inQ,
                        const SparseIndices &p1q2, const SparseIndices &p2q1) {
   SparseIndices p1q1(3 * p1q2.size() + 3 * p2q1.size());
-  thrust::for_each_n(thrust::device, zip(countAt(0), p1q2.begin(0), p1q2.begin(1)),
+  auto policy = autoPolicy(p1q2.size());
+  for_each_n(policy, zip(countAt(0), p1q2.begin(0), p1q2.begin(1)),
                      p1q2.size(),
                      CopyFaceEdges({p1q1.ptrDpq(), inQ.halfedge_.cptrD()}));
 
   p1q1.SwapPQ();
-  thrust::for_each_n(thrust::device, zip(countAt(p1q2.size()), p2q1.begin(1), p2q1.begin(0)),
+  for_each_n(policy, zip(countAt(p1q2.size()), p2q1.begin(1), p2q1.begin(0)),
                      p2q1.size(),
                      CopyFaceEdges({p1q1.ptrDpq(), inP.halfedge_.cptrD()}));
   p1q1.SwapPQ();
@@ -245,8 +247,8 @@ std::tuple<VecDH<int>, VecDH<glm::vec4>> Shadow11(SparseIndices &p1q1,
   VecDH<int> s11(p1q1.size());
   VecDH<glm::vec4> xyzz11(p1q1.size());
 
-  thrust::for_each_n(
-      thrust::device, zip(xyzz11.begin(), s11.begin(), p1q1.begin(0), p1q1.begin(1)),
+  for_each_n(
+      autoPolicy(p1q1.size()), zip(xyzz11.begin(), s11.begin(), p1q1.begin(0), p1q1.begin(1)),
       p1q1.size(),
       Kernel11({inP.vertPos_.cptrD(), inQ.vertPos_.cptrD(),
                 inP.halfedge_.cptrD(), inQ.halfedge_.cptrD(), expandP,
@@ -342,8 +344,8 @@ std::tuple<VecDH<int>, VecDH<float>> Shadow02(const Manifold::Impl &inP,
 
   auto vertNormalP =
       forward ? inP.vertNormal_.cptrD() : inQ.vertNormal_.cptrD();
-  thrust::for_each_n(
-      thrust::device, zip(s02.begin(), z02.begin(), p0q2.begin(!forward),
+  for_each_n(
+      autoPolicy(p0q2.size()), zip(s02.begin(), z02.begin(), p0q2.begin(!forward),
           p0q2.begin(forward)),
       p0q2.size(),
       Kernel02({inP.vertPos_.cptrD(), inQ.halfedge_.cptrD(),
@@ -452,8 +454,8 @@ std::tuple<VecDH<int>, VecDH<glm::vec3>> Intersect12(
   VecDH<int> x12(p1q2.size());
   VecDH<glm::vec3> v12(p1q2.size());
 
-  thrust::for_each_n(
-      thrust::device, zip(x12.begin(), v12.begin(), p1q2.begin(!forward),
+  for_each_n(
+      autoPolicy(p1q2.size()), zip(x12.begin(), v12.begin(), p1q2.begin(!forward),
           p1q2.begin(forward)),
       p1q2.size(),
       Kernel12({p0q2.ptrDpq(), s02.ptrD(), z02.cptrD(), p0q2.size(),
@@ -471,19 +473,21 @@ VecDH<int> Winding03(const Manifold::Impl &inP, SparseIndices &p0q2,
   // verts that are not shadowed (not in p0q2) have winding number zero.
   VecDH<int> w03(inP.NumVert(), 0);
 
-  if (!thrust::is_sorted(thrust::device, p0q2.begin(reverse), p0q2.end(reverse)))
-    thrust::sort_by_key(thrust::device, p0q2.begin(reverse), p0q2.end(reverse), s02.begin());
+  auto policy = autoPolicy(p0q2.size());
+  if (!is_sorted(policy, p0q2.begin(reverse), p0q2.end(reverse)))
+    sort_by_key(policy, p0q2.begin(reverse), p0q2.end(reverse), s02.begin());
   VecDH<int> w03val(w03.size());
   VecDH<int> w03vert(w03.size());
   // sum known s02 values into w03 (winding number)
   auto endPair =
-      thrust::reduce_by_key(thrust::device, p0q2.begin(reverse), p0q2.end(reverse),
+      reduce_by_key<thrust::pair<decltype(w03val.begin()), decltype(w03val.begin())>>
+      (policy, p0q2.begin(reverse), p0q2.end(reverse),
                             s02.begin(), w03vert.begin(), w03val.begin());
-  thrust::scatter(thrust::device, w03val.begin(), endPair.second, w03vert.begin(),
+  scatter(autoPolicy(endPair.second - w03val.begin()), w03val.begin(), endPair.second, w03vert.begin(),
                   w03.begin());
 
   if (reverse)
-    thrust::transform(thrust::device, w03.begin(), w03.end(), w03.begin(),
+    transform(autoPolicy(w03.size()), w03.begin(), w03.end(), w03.begin(),
                       thrust::negate<int>());
   return w03;
 };

manifold/src/boolean_result.cpp

@@ -17,6 +17,7 @@
 
 #include "boolean3.h"
 #include "polygon.h"
+#include "par.h"
 
 // TODO: make this runtime configurable for quicker debug
 constexpr bool kVerbose = false;
@@ -79,49 +80,51 @@ std::tuple<VecDH<int>, VecDH<int>> SizeOutput(
     const VecDH<int> &i03, const VecDH<int> &i30, const VecDH<int> &i12,
     const VecDH<int> &i21, const SparseIndices &p1q2, const SparseIndices &p2q1,
     bool invertQ) {
-  VecDH<int> sidesPerFacePQ(inP.NumTri() + inQ.NumTri());
+  VecDH<int> sidesPerFacePQ(inP.NumTri() + inQ.NumTri(), 0);
   auto sidesPerFaceP = sidesPerFacePQ.ptrD();
   auto sidesPerFaceQ = sidesPerFacePQ.ptrD() + inP.NumTri();
 
-  thrust::for_each(thrust::device, inP.halfedge_.begin(), inP.halfedge_.end(),
+  auto policy = autoPolicy(std::max(inP.halfedge_.size(), inQ.halfedge_.size()));
+
+  for_each(policy, inP.halfedge_.begin(), inP.halfedge_.end(),
                    CountVerts({sidesPerFaceP, i03.cptrD()}));
-  thrust::for_each(thrust::device, inQ.halfedge_.begin(), inQ.halfedge_.end(),
+  for_each(policy, inQ.halfedge_.begin(), inQ.halfedge_.end(),
                    CountVerts({sidesPerFaceQ, i30.cptrD()}));
-  thrust::for_each_n(
-      thrust::device, zip(p1q2.begin(0), p1q2.begin(1), i12.begin()), i12.size(),
+  for_each_n(
+      policy, zip(p1q2.begin(0), p1q2.begin(1), i12.begin()), i12.size(),
       CountNewVerts({sidesPerFaceP, sidesPerFaceQ, inP.halfedge_.cptrD()}));
-  thrust::for_each_n(
-      thrust::device, zip(p2q1.begin(1), p2q1.begin(0), i21.begin()), i21.size(),
+  for_each_n(
+      policy, zip(p2q1.begin(1), p2q1.begin(0), i21.begin()), i21.size(),
       CountNewVerts({sidesPerFaceQ, sidesPerFaceP, inQ.halfedge_.cptrD()}));
 
-  VecDH<int> facePQ2R(inP.NumTri() + inQ.NumTri() + 1);
+  VecDH<int> facePQ2R(inP.NumTri() + inQ.NumTri() + 1, 0);
   auto keepFace =
       thrust::make_transform_iterator(sidesPerFacePQ.begin(), NotZero());
-  thrust::inclusive_scan(thrust::device, keepFace, keepFace + sidesPerFacePQ.size(),
+  inclusive_scan(policy, keepFace, keepFace + sidesPerFacePQ.size(),
                          facePQ2R.begin() + 1);
   int numFaceR = facePQ2R.back();
   facePQ2R.resize(inP.NumTri() + inQ.NumTri());
 
   outR.faceNormal_.resize(numFaceR);
-  auto next = thrust::copy_if(thrust::device, inP.faceNormal_.begin(), inP.faceNormal_.end(),
+  auto next = copy_if<decltype(outR.faceNormal_.begin())>(policy, inP.faceNormal_.begin(), inP.faceNormal_.end(),
                               keepFace, outR.faceNormal_.begin(),
                               thrust::identity<bool>());
   if (invertQ) {
     auto start = thrust::make_transform_iterator(inQ.faceNormal_.begin(),
                                                  thrust::negate<glm::vec3>());
     auto end = thrust::make_transform_iterator(inQ.faceNormal_.end(),
                                                thrust::negate<glm::vec3>());
-    thrust::copy_if(thrust::device, start, end, keepFace + inP.NumTri(), next,
+    copy_if<decltype(inQ.faceNormal_.begin())>(policy, start, end, keepFace + inP.NumTri(), next,
                     thrust::identity<bool>());
   } else {
-    thrust::copy_if(thrust::device, inQ.faceNormal_.begin(), inQ.faceNormal_.end(),
+    copy_if<decltype(inQ.faceNormal_.begin())>(policy, inQ.faceNormal_.begin(), inQ.faceNormal_.end(),
                     keepFace + inP.NumTri(), next, thrust::identity<bool>());
   }
 
   auto newEnd =
-      thrust::remove(thrust::device, sidesPerFacePQ.begin(), sidesPerFacePQ.end(), 0);
-  VecDH<int> faceEdge(newEnd - sidesPerFacePQ.begin() + 1);
-  thrust::inclusive_scan(thrust::device, sidesPerFacePQ.begin(), newEnd,
+      remove<decltype(sidesPerFacePQ.begin())>(policy, sidesPerFacePQ.begin(), sidesPerFacePQ.end(), 0);
+  VecDH<int> faceEdge(newEnd - sidesPerFacePQ.begin() + 1, 0);
+  inclusive_scan(policy, sidesPerFacePQ.begin(), newEnd,
                          faceEdge.begin() + 1);
   outR.halfedge_.resize(faceEdge.back());
 
@@ -405,8 +408,8 @@ void AppendWholeEdges(Manifold::Impl &outR, VecDH<int> &facePtrR,
                       const VecDH<char> wholeHalfedgeP, const VecDH<int> &i03,
                       const VecDH<int> &vP2R, const int *faceP2R,
                       bool forward) {
-  thrust::for_each_n(
-      thrust::device, zip(wholeHalfedgeP.begin(), inP.halfedge_.begin(), countAt(0)),
+  for_each_n(
+      autoPolicy(inP.halfedge_.size()), zip(wholeHalfedgeP.begin(), inP.halfedge_.begin(), countAt(0)),
       inP.halfedge_.size(),
       DuplicateHalfedges({outR.halfedge_.ptrD(), halfedgeRef.ptrD(),
                           facePtrR.ptrD(), inP.halfedge_.cptrD(), i03.cptrD(),
@@ -481,8 +484,8 @@ std::pair<VecDH<BaryRef>, VecDH<int>> CalculateMeshRelation(
   VecDH<BaryRef> faceRef(numFaceR);
   VecDH<int> halfedgeBary(halfedgeRef.size());
   VecDH<int> idx(1, 0);
-  thrust::for_each_n(
-      thrust::device, zip(halfedgeBary.begin(), halfedgeRef.begin(),
+  for_each_n(
+      autoPolicy(halfedgeRef.size()), zip(halfedgeBary.begin(), halfedgeRef.begin(),
           outR.halfedge_.cbegin()),
       halfedgeRef.size(),
       CreateBarycentric(
@@ -550,33 +553,35 @@ Manifold::Impl Boolean3::Result(Manifold::OpType op) const {
   VecDH<int> i21(x21_.size());
   VecDH<int> i03(w03_.size());
   VecDH<int> i30(w30_.size());
-  thrust::transform(thrust::device, x12_.begin(), x12_.end(), i12.begin(), c3 * _1);
-  thrust::transform(thrust::device, x21_.begin(), x21_.end(), i21.begin(), c3 * _1);
-  thrust::transform(thrust::device, w03_.begin(), w03_.end(), i03.begin(), c1 + c3 * _1);
-  thrust::transform(thrust::device, w30_.begin(), w30_.end(), i30.begin(), c2 + c3 * _1);
+  auto policy = autoPolicy(std::max(std::max(x12_.size(), x21_.size()), std::max(w03_.size(), w30_.size())));
+
+  transform(policy, x12_.begin(), x12_.end(), i12.begin(), c3 * _1);
+  transform(policy, x21_.begin(), x21_.end(), i21.begin(), c3 * _1);
+  transform(policy, w03_.begin(), w03_.end(), i03.begin(), c1 + c3 * _1);
+  transform(policy, w30_.begin(), w30_.end(), i30.begin(), c2 + c3 * _1);
 
   VecDH<int> vP2R(inP_.NumVert());
-  thrust::exclusive_scan(thrust::device, i03.begin(), i03.end(), vP2R.begin(), 0, AbsSum());
+  exclusive_scan(policy, i03.begin(), i03.end(), vP2R.begin(), 0, AbsSum());
   int numVertR = AbsSum()(vP2R.back(), i03.back());
   const int nPv = numVertR;
 
   VecDH<int> vQ2R(inQ_.NumVert());
-  thrust::exclusive_scan(thrust::device, i30.begin(), i30.end(), vQ2R.begin(), numVertR,
+  exclusive_scan(policy, i30.begin(), i30.end(), vQ2R.begin(), numVertR,
                          AbsSum());
   numVertR = AbsSum()(vQ2R.back(), i30.back());
   const int nQv = numVertR - nPv;
 
   VecDH<int> v12R(v12_.size());
   if (v12_.size() > 0) {
-    thrust::exclusive_scan(thrust::device, i12.begin(), i12.end(), v12R.begin(), numVertR,
+    exclusive_scan(policy, i12.begin(), i12.end(), v12R.begin(), numVertR,
                            AbsSum());
     numVertR = AbsSum()(v12R.back(), i12.back());
   }
   const int n12 = numVertR - nPv - nQv;
 
   VecDH<int> v21R(v21_.size());
   if (v21_.size() > 0) {
-    thrust::exclusive_scan(thrust::device, i21.begin(), i21.end(), v21R.begin(), numVertR,
+    exclusive_scan(policy, i21.begin(), i21.end(), v21R.begin(), numVertR,
                            AbsSum());
     numVertR = AbsSum()(v21R.back(), i21.back());
   }
@@ -592,14 +597,14 @@ Manifold::Impl Boolean3::Result(Manifold::OpType op) const {
   outR.vertPos_.resize(numVertR);
   // Add vertices, duplicating for inclusion numbers not in [-1, 1].
   // Retained vertices from P and Q:
-  thrust::for_each_n(thrust::device, zip(i03.begin(), vP2R.begin(), inP_.vertPos_.begin()),
+  for_each_n(policy, zip(i03.begin(), vP2R.begin(), inP_.vertPos_.begin()),
                      inP_.NumVert(), DuplicateVerts({outR.vertPos_.ptrD()}));
-  thrust::for_each_n(thrust::device, zip(i30.begin(), vQ2R.begin(), inQ_.vertPos_.begin()),
+  for_each_n(policy, zip(i30.begin(), vQ2R.begin(), inQ_.vertPos_.begin()),
                      inQ_.NumVert(), DuplicateVerts({outR.vertPos_.ptrD()}));
   // New vertices created from intersections:
-  thrust::for_each_n(thrust::device, zip(i12.begin(), v12R.begin(), v12_.begin()),
+  for_each_n(policy, zip(i12.begin(), v12R.begin(), v12_.begin()),
                      i12.size(), DuplicateVerts({outR.vertPos_.ptrD()}));
-  thrust::for_each_n(thrust::device, zip(i21.begin(), v21R.begin(), v21_.begin()),
+  for_each_n(policy, zip(i21.begin(), v21R.begin(), v21_.begin()),
                      i21.size(), DuplicateVerts({outR.vertPos_.ptrD()}));
 
   if (kVerbose) {