diff --git a/.github/workflows/manifold.yml b/.github/workflows/manifold.yml index 010a0534b..3daeb520b 100644 --- a/.github/workflows/manifold.yml +++ b/.github/workflows/manifold.yml @@ -221,7 +221,7 @@ jobs: pip install trimesh pytest - name: Install TBB if: matrix.parallel_backend == 'TBB' - run: brew install tbb + run: brew install tbb onedpl - uses: actions/checkout@v4 with: submodules: recursive diff --git a/.vscode/settings.json b/.vscode/settings.json index 378812de1..bb7698552 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -113,7 +113,8 @@ "format": "cpp", "numbers": "cpp", "ranges": "cpp", - "stop_token": "cpp" + "stop_token": "cpp", + "execution": "cpp" }, "C_Cpp.clang_format_fallbackStyle": "google", "editor.formatOnSave": true, diff --git a/README.md b/README.md index 37897d705..c30ce70a6 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,6 @@ This is a modern C++ library that Github's CI verifies builds and runs on a vari System Dependencies (note that we will automatically download the dependency if there is no such package on the system): - [`GLM`](https://github.com/g-truc/glm/): A compact header-only vector library. -- [`Thrust`](https://github.com/NVIDIA/thrust): NVIDIA's parallel algorithms library (basically a superset of C++17 std::parallel_algorithms) - [`tbb`](https://github.com/oneapi-src/oneTBB/): Intel's thread building blocks library. (only when `MANIFOLD_PAR=TBB` is enabled) - [`gtest`](https://github.com/google/googletest/): Google test library (only when test is enabled, i.e. `MANIFOLD_TEST=ON`) @@ -48,7 +47,7 @@ This library is fast with guaranteed manifold output. As such you need manifold The most significant contribution here is a guaranteed-manifold [mesh Boolean](https://github.com/elalish/manifold/wiki/Manifold-Library#mesh-boolean) algorithm, which I believe is the first of its kind. If you know of another, please open a discussion - a mesh Boolean algorithm robust to edge cases has been an open problem for many years. Likewise, if the Boolean here ever fails you, please submit an issue! This Boolean forms the basis of a CAD kernel, as it allows simple shapes to be combined into more complex ones. -To aid in speed, this library makes extensive use of parallelization, generally through Nvidia's Thrust library. You can switch between the TBB, and serial C++ backends by setting a CMake flag. Not everything is so parallelizable, for instance a [polygon triangulation](https://github.com/elalish/manifold/wiki/Manifold-Library#polygon-triangulation) algorithm is included which is serial. Even if compiled with parallel backend, the code will still fall back to the serial version of the algorithms if the problem size is small. The WASM build is serial-only for now, but still fast. +To aid in speed, this library makes extensive use of parallelization, generally through PSTL. You can switch between the TBB, and serial C++ backends by setting a CMake flag. Not everything is so parallelizable, for instance a [polygon triangulation](https://github.com/elalish/manifold/wiki/Manifold-Library#polygon-triangulation) algorithm is included which is serial. Even if compiled with parallel backend, the code will still fall back to the serial version of the algorithms if the problem size is small. The WASM build is serial-only for now, but still fast. > Note: OMP and CUDA backends are now removed @@ -83,7 +82,6 @@ CMake flags (usage e.g. `-DMANIFOLD_DEBUG=ON`): Offline building: - `FETCHCONTENT_SOURCE_DIR_GLM`: path to glm source. - `FETCHCONTENT_SOURCE_DIR_GOOGLETEST`: path to googletest source. -- `FETCHCONTENT_SOURCE_DIR_THRUST`: path to NVIDIA thrust source. The build instructions used by our CI are in [manifold.yml](https://github.com/elalish/manifold/blob/master/.github/workflows/manifold.yml), which is a good source to check if something goes wrong and for instructions specific to other platforms, like Windows. diff --git a/flake.lock b/flake.lock index de866f1f1..a9dad5445 100644 --- a/flake.lock +++ b/flake.lock @@ -71,8 +71,7 @@ "clipper2-src": "clipper2-src", "flake-utils": "flake-utils", "gtest-src": "gtest-src", - "nixpkgs": "nixpkgs", - "thrust-src": "thrust-src" + "nixpkgs": "nixpkgs" } }, "systems": { @@ -89,24 +88,6 @@ "repo": "default", "type": "github" } - }, - "thrust-src": { - "flake": false, - "locked": { - "lastModified": 1696873248, - "narHash": "sha256-mUhMXGPbO2t83EvI8YNDRLngvLiwfzs4EgRrcxKfhHs=", - "ref": "refs/heads/main", - "rev": "756c5afc0750f1413da05bd2b6505180e84c53d4", - "revCount": 4758, - "submodules": true, - "type": "git", - "url": "https://github.com/NVIDIA/thrust.git" - }, - "original": { - "submodules": true, - "type": "git", - "url": "https://github.com/NVIDIA/thrust.git" - } } }, "root": "root", diff --git a/flake.nix b/flake.nix index c07cf3b77..05677219a 100644 --- a/flake.nix +++ b/flake.nix @@ -5,15 +5,11 @@ url = "github:google/googletest/v1.14.0"; flake = false; }; - inputs.thrust-src = { - url = "git+https://github.com/NVIDIA/thrust.git?submodules=1"; - flake = false; - }; inputs.clipper2-src = { url = "github:AngusJohnson/Clipper2"; flake = false; }; - outputs = { self, nixpkgs, flake-utils, gtest-src, thrust-src, clipper2-src }: + outputs = { self, nixpkgs, flake-utils, gtest-src, clipper2-src }: flake-utils.lib.eachDefaultSystem (system: let @@ -52,7 +48,6 @@ "-DMANIFOLD_PYBIND=ON" "-DMANIFOLD_CBIND=ON" "-DBUILD_SHARED_LIBS=ON" - "-DFETCHCONTENT_SOURCE_DIR_THRUST=${thrust-src}" "-DMANIFOLD_PAR=${pkgs.lib.strings.toUpper parallel-backend}" ]; prePatch = '' @@ -100,7 +95,6 @@ emcmake cmake -DCMAKE_BUILD_TYPE=Release \ -DFETCHCONTENT_SOURCE_DIR_GLM=${pkgs.glm.src} \ -DFETCHCONTENT_SOURCE_DIR_GOOGLETEST=${gtest-src} \ - -DFETCHCONTENT_SOURCE_DIR_THRUST=${thrust-src} \ -DFETCHCONTENT_SOURCE_DIR_CLIPPER2=../clipper2 .. ''; buildPhase = '' @@ -139,7 +133,6 @@ pathspec pkg-config ]; - SKBUILD_CMAKE_DEFINE = "FETCHCONTENT_SOURCE_DIR_THRUST=${thrust-src}"; checkInputs = [ trimesh pytest diff --git a/manifoldDeps.cmake b/manifoldDeps.cmake index 9dcbaccaf..e0d844665 100644 --- a/manifoldDeps.cmake +++ b/manifoldDeps.cmake @@ -1,9 +1,12 @@ include(FetchContent) include(GNUInstallDirs) find_package(PkgConfig QUIET) -find_package(Clipper2) +find_package(Clipper2 QUIET) if(MANIFOLD_PAR STREQUAL "TBB") find_package(TBB QUIET) + if(APPLE) + find_package(oneDPL QUIET) + endif() endif() if (PKG_CONFIG_FOUND) if (NOT Clipper2_FOUND) diff --git a/src/collider/include/collider.h b/src/collider/include/collider.h index bf5235445..8ffed384c 100644 --- a/src/collider/include/collider.h +++ b/src/collider/include/collider.h @@ -36,7 +36,7 @@ class Collider { Vec nodeBBox_; Vec nodeParent_; // even nodes are leaves, odd nodes are internal, root is 1 - Vec> internalChildren_; + Vec> internalChildren_; size_t NumInternal() const { return internalChildren_.size(); }; size_t NumLeaves() const { diff --git a/src/collider/src/collider.cpp b/src/collider/src/collider.cpp index 02c1c35bb..9a9248363 100644 --- a/src/collider/src/collider.cpp +++ b/src/collider/src/collider.cpp @@ -69,7 +69,7 @@ int Leaf2Node(int leaf) { return leaf * 2; } struct CreateRadixTree { VecView nodeParent_; - VecView> internalChildren_; + VecView> internalChildren_; const VecView leafMorton_; int PrefixLength(uint32_t a, uint32_t b) const { @@ -137,7 +137,7 @@ struct CreateRadixTree { int first = internal; // Find the range of objects with a common prefix int last = RangeEnd(first); - if (first > last) thrust::swap(first, last); + if (first > last) std::swap(first, last); // Determine where the next-highest difference occurs int split = FindSplit(first, last); int child1 = split == first ? Leaf2Node(split) : Internal2Node(split); @@ -156,7 +156,7 @@ template struct FindCollisions { VecView queries; VecView nodeBBox_; - VecView> internalChildren_; + VecView> internalChildren_; Recorder recorder; int RecordCollision(int node, const int queryIdx) { @@ -244,7 +244,7 @@ struct BuildInternalBoxes { VecView nodeBBox_; VecView counter_; const VecView nodeParent_; - const VecView> internalChildren_; + const VecView> internalChildren_; void operator()(int leaf) { int node = Leaf2Node(leaf); @@ -288,7 +288,7 @@ Collider::Collider(const VecView& leafBB, // assign and allocate members nodeBBox_.resize(num_nodes); nodeParent_.resize(num_nodes, -1); - internalChildren_.resize(leafBB.size() - 1, thrust::make_pair(-1, -1)); + internalChildren_.resize(leafBB.size() - 1, std::make_pair(-1, -1)); // organize tree for_each_n(autoPolicy(NumInternal()), countAt(0), NumInternal(), CreateRadixTree({nodeParent_, internalChildren_, leafMorton})); diff --git a/src/manifold/src/boolean3.cpp b/src/manifold/src/boolean3.cpp index e80cd0ec6..eb3c03afa 100644 --- a/src/manifold/src/boolean3.cpp +++ b/src/manifold/src/boolean3.cpp @@ -105,7 +105,7 @@ inline bool Shadows(float p, float q, float dir) { return p == q ? dir < 0 : p < q; } -inline thrust::pair Shadow01( +inline std::pair Shadow01( const int p0, const int q1, VecView vertPosP, VecView vertPosQ, VecView halfedgeQ, const float expandP, VecView normalP, const bool reverse) { @@ -133,7 +133,7 @@ inline thrust::pair Shadow01( if (!Shadows(vertPosP[p0].y, yz01[0], expandP * normalP[p0].y)) s01 = 0; } } - return thrust::make_pair(s01, yz01); + return std::make_pair(s01, yz01); } // https://github.com/scandum/binary_search/blob/master/README.md @@ -406,7 +406,7 @@ struct Kernel12 { if (k < 2 && (k == 0 || (s != 0) != shadows)) { shadows = s != 0; xzyLR0[k] = vertPosP[vert]; - thrust::swap(xzyLR0[k].y, xzyLR0[k].z); + std::swap(xzyLR0[k].y, xzyLR0[k].z); xzyLR1[k] = xzyLR0[k]; xzyLR1[k][1] = z02[idx]; k++; @@ -434,7 +434,7 @@ struct Kernel12 { xzyLR0[k][2] = xyzz.y; xzyLR1[k] = xzyLR0[k]; xzyLR1[k][1] = xyzz.w; - if (!forward) thrust::swap(xzyLR0[k][1], xzyLR1[k][1]); + if (!forward) std::swap(xzyLR0[k][1], xzyLR1[k][1]); k++; } } diff --git a/src/manifold/src/boolean_result.cpp b/src/manifold/src/boolean_result.cpp index 3f6760538..8eea76201 100644 --- a/src/manifold/src/boolean_result.cpp +++ b/src/manifold/src/boolean_result.cpp @@ -14,6 +14,7 @@ #include #include +#include #include #if MANIFOLD_PAR == 'T' && __has_include() @@ -30,10 +31,8 @@ using concurrent_map = std::map; #endif #include "boolean3.h" #include "par.h" -#include "polygon.h" using namespace manifold; -using namespace thrust::placeholders; template <> struct std::hash> { @@ -46,7 +45,7 @@ namespace { constexpr int kParallelThreshold = 128; -struct AbsSum : public thrust::binary_function { +struct AbsSum { int operator()(int a, int b) { return abs(a) + abs(b); } }; @@ -93,10 +92,6 @@ struct CountNewVerts { } }; -struct NotZero : public thrust::unary_function { - int operator()(int x) const { return x > 0 ? 1 : 0; } -}; - std::tuple, Vec> SizeOutput( Manifold::Impl &outR, const Manifold::Impl &inP, const Manifold::Impl &inQ, const Vec &i03, const Vec &i30, const Vec &i12, @@ -122,7 +117,8 @@ std::tuple, Vec> SizeOutput( {sidesPerFaceQ, sidesPerFaceP, i21, p2q1, inQ.halfedge_})); Vec facePQ2R(inP.NumTri() + inQ.NumTri() + 1, 0); - auto keepFace = TransformIterator(sidesPerFacePQ.begin(), NotZero()); + auto keepFace = TransformIterator(sidesPerFacePQ.begin(), + [](int x) { return x > 0 ? 1 : 0; }); inclusive_scan(autoPolicy(sidesPerFacePQ.size()), keepFace, keepFace + sidesPerFacePQ.size(), facePQ2R.begin() + 1); @@ -130,23 +126,39 @@ std::tuple, Vec> SizeOutput( facePQ2R.resize(inP.NumTri() + inQ.NumTri()); outR.faceNormal_.resize(numFaceR); - auto next = copy_if( - autoPolicy(inP.faceNormal_.size()), inP.faceNormal_.begin(), - inP.faceNormal_.end(), keepFace, outR.faceNormal_.begin(), - thrust::identity()); + + Vec tmpBuffer(outR.faceNormal_.size()); + auto faceIds = TransformIterator(countAt(0_z), [&sidesPerFacePQ](size_t i) { + if (sidesPerFacePQ[i] > 0) return i; + return std::numeric_limits::max(); + }); + + auto next = copy_if( + autoPolicy(inP.faceNormal_.size()), faceIds, + faceIds + inP.faceNormal_.size(), tmpBuffer.begin(), + [](size_t v) { return v != std::numeric_limits::max(); }); + + gather(autoPolicy(inP.faceNormal_.size()), tmpBuffer.begin(), next, + inP.faceNormal_.begin(), outR.faceNormal_.begin()); + + auto faceIdsQ = + TransformIterator(countAt(0_z), [&sidesPerFacePQ, &inP](size_t i) { + if (sidesPerFacePQ[i + inP.faceNormal_.size()] > 0) return i; + return std::numeric_limits::max(); + }); + auto end = copy_if( + autoPolicy(inQ.faceNormal_.size()), faceIdsQ, + faceIdsQ + inQ.faceNormal_.size(), next, + [](size_t v) { return v != std::numeric_limits::max(); }); + if (invertQ) { - auto start = - TransformIterator(inQ.faceNormal_.begin(), thrust::negate()); - auto end = - TransformIterator(inQ.faceNormal_.end(), thrust::negate()); - copy_if( - autoPolicy(inQ.faceNormal_.size()), start, end, keepFace + inP.NumTri(), - next, thrust::identity()); + gather(autoPolicy(inQ.faceNormal_.size()), next, end, + TransformIterator(inQ.faceNormal_.begin(), Negate()), + outR.faceNormal_.begin() + std::distance(tmpBuffer.begin(), next)); } else { - copy_if( - autoPolicy(inQ.faceNormal_.size()), inQ.faceNormal_.begin(), - inQ.faceNormal_.end(), keepFace + inP.NumTri(), next, - thrust::identity()); + gather(autoPolicy(inQ.faceNormal_.size()), next, end, + inQ.faceNormal_.begin(), + outR.faceNormal_.begin() + std::distance(tmpBuffer.begin(), next)); } auto newEnd = remove( @@ -662,13 +674,13 @@ Manifold::Impl Boolean3::Result(OpType op) const { Vec i30(w30_.size()); transform(autoPolicy(x12_.size()), x12_.begin(), x12_.end(), i12.begin(), - c3 * _1); + [c3](int v) { return c3 * v; }); transform(autoPolicy(x21_.size()), x21_.begin(), x21_.end(), i21.begin(), - c3 * _1); + [c3](int v) { return c3 * v; }); transform(autoPolicy(w03_.size()), w03_.begin(), w03_.end(), i03.begin(), - c1 + c3 * _1); + [c1, c3](int v) { return c1 + c3 * v; }); transform(autoPolicy(w30_.size()), w30_.begin(), w30_.end(), i30.begin(), - c2 + c3 * _1); + [c2, c3](int v) { return c2 + c3 * v; }); Vec vP2R(inP_.NumVert()); exclusive_scan(autoPolicy(i03.size()), i03.begin(), i03.end(), vP2R.begin(), diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp index e1e891b4d..a2975c8f6 100644 --- a/src/manifold/src/impl.cpp +++ b/src/manifold/src/impl.cpp @@ -129,8 +129,8 @@ struct UpdateMeshID { }; struct CoplanarEdge { - VecView> face2face; - VecView> vert2vert; + VecView> face2face; + VecView> vert2vert; VecView triArea; VecView halfedge; VecView vertPos; @@ -164,7 +164,7 @@ struct CoplanarEdge { } } if (propEqual) { - vert2vert[edgeIdx] = thrust::make_pair(prop0, prop1); + vert2vert[edgeIdx] = std::make_pair(prop0, prop1); } } @@ -221,7 +221,7 @@ struct CoplanarEdge { } } - face2face[edgeIdx] = thrust::make_pair(edge.face, pair.face); + face2face[edgeIdx] = std::make_pair(edge.face, pair.face); } }; @@ -257,7 +257,7 @@ struct CheckCoplanarity { }; int GetLabels(std::vector& components, - const Vec>& edges, int numNodes) { + const Vec>& edges, int numNodes) { UnionFind<> uf(numNodes); for (auto edge : edges) { if (edge.first == -1 || edge.second == -1) continue; @@ -268,7 +268,7 @@ int GetLabels(std::vector& components, } void DedupePropVerts(manifold::Vec& triProp, - const Vec>& vert2vert) { + const Vec>& vert2vert) { ZoneScoped; std::vector vertLabels; const int numLabels = GetLabels(vertLabels, vert2vert, vert2vert.size()); @@ -539,11 +539,21 @@ void Manifold::Impl::RemoveUnreferencedVerts(Vec& triVerts) { MarkVerts({vertOld2New.view(1)})); const Vec oldVertPos = vertPos_; - vertPos_.resize(copy_if( - policy, oldVertPos.cbegin(), oldVertPos.cend(), - vertOld2New.cbegin() + 1, vertPos_.begin(), - thrust::identity()) - - vertPos_.begin()); + + Vec tmpBuffer(oldVertPos.size()); + auto vertIdIter = TransformIterator(countAt(0_z), [&vertOld2New](size_t i) { + if (vertOld2New[i + 1] > 0) return i; + return std::numeric_limits::max(); + }); + + auto next = copy_if( + autoPolicy(tmpBuffer.size()), vertIdIter, vertIdIter + tmpBuffer.size(), + tmpBuffer.begin(), + [](size_t v) { return v != std::numeric_limits::max(); }); + gather(autoPolicy(tmpBuffer.size()), tmpBuffer.begin(), next, + oldVertPos.begin(), vertPos_.begin()); + + vertPos_.resize(std::distance(tmpBuffer.begin(), next)); inclusive_scan(policy, vertOld2New.begin() + 1, vertOld2New.end(), vertOld2New.begin() + 1); @@ -572,8 +582,8 @@ void Manifold::Impl::CreateFaces(const std::vector& propertyTolerance) { propertyTolerance.empty() ? Vec(meshRelation_.numProp, kTolerance) : propertyTolerance; - Vec> face2face(halfedge_.size(), {-1, -1}); - Vec> vert2vert(halfedge_.size(), {-1, -1}); + Vec> face2face(halfedge_.size(), {-1, -1}); + Vec> vert2vert(halfedge_.size(), {-1, -1}); Vec triArea(NumTri()); for_each_n(autoPolicy(halfedge_.size()), countAt(0), halfedge_.size(), CoplanarEdge({face2face, vert2vert, triArea, halfedge_, vertPos_, @@ -680,7 +690,7 @@ void Manifold::Impl::MarkFailure(Error status) { void Manifold::Impl::Warp(std::function warpFunc) { WarpBatch([&warpFunc](VecView vecs) { - thrust::for_each(thrust::host, vecs.begin(), vecs.end(), warpFunc); + for_each(ExecutionPolicy::Seq, vecs.begin(), vecs.end(), warpFunc); }); } diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index e24aef2e2..d7d44f1f8 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -21,11 +21,9 @@ #include "csg_tree.h" #include "impl.h" #include "par.h" -#include "tri_dist.h" namespace { using namespace manifold; -using namespace thrust::placeholders; ExecutionParams manifoldParams; @@ -612,12 +610,11 @@ Manifold Manifold::SetProperties( } else { pImpl->meshRelation_.properties = Vec(numProp * NumPropVert(), 0); } - thrust::for_each_n( - thrust::host, countAt(0), NumTri(), - UpdateProperties({pImpl->meshRelation_.properties.data(), numProp, - oldProperties.data(), oldNumProp, - pImpl->vertPos_.data(), triProperties.data(), - pImpl->halfedge_.data(), propFunc})); + for_each_n(ExecutionPolicy::Seq, countAt(0), NumTri(), + UpdateProperties({pImpl->meshRelation_.properties.data(), + numProp, oldProperties.data(), oldNumProp, + pImpl->vertPos_.data(), triProperties.data(), + pImpl->halfedge_.data(), propFunc})); } pImpl->meshRelation_.numProp = numProp; diff --git a/src/manifold/src/mesh_fixes.h b/src/manifold/src/mesh_fixes.h index 00a2782df..461bef93b 100644 --- a/src/manifold/src/mesh_fixes.h +++ b/src/manifold/src/mesh_fixes.h @@ -39,11 +39,10 @@ struct FlipTris { VecView halfedge; void operator()(const int tri) { - thrust::swap(halfedge[3 * tri], halfedge[3 * tri + 2]); + std::swap(halfedge[3 * tri], halfedge[3 * tri + 2]); for (const int i : {0, 1, 2}) { - thrust::swap(halfedge[3 * tri + i].startVert, - halfedge[3 * tri + i].endVert); + std::swap(halfedge[3 * tri + i].startVert, halfedge[3 * tri + i].endVert); halfedge[3 * tri + i].pairedHalfedge = FlipHalfedge(halfedge[3 * tri + i].pairedHalfedge); } diff --git a/src/manifold/src/properties.cpp b/src/manifold/src/properties.cpp index 6bc62fe5e..082d84c7b 100644 --- a/src/manifold/src/properties.cpp +++ b/src/manifold/src/properties.cpp @@ -26,7 +26,7 @@ struct FaceAreaVolume { VecView vertPos; const float precision; - thrust::pair operator()(int face) { + std::pair operator()(int face) { float perimeter = 0; glm::vec3 edge[3]; for (int i : {0, 1, 2}) { @@ -40,56 +40,7 @@ struct FaceAreaVolume { float area = glm::length(crossP); float volume = glm::dot(crossP, vertPos[halfedges[3 * face].startVert]); - return thrust::make_pair(area / 2.0f, volume / 6.0f); - } -}; - -struct PosMin - : public thrust::binary_function { - glm::vec3 operator()(glm::vec3 a, glm::vec3 b) { - if (isnan(a.x)) return b; - if (isnan(b.x)) return a; - return glm::min(a, b); - } -}; - -struct PosMax - : public thrust::binary_function { - glm::vec3 operator()(glm::vec3 a, glm::vec3 b) { - if (isnan(a.x)) return b; - if (isnan(b.x)) return a; - return glm::max(a, b); - } -}; - -struct FiniteVert { - bool operator()(glm::vec3 v) { return glm::all(glm::isfinite(v)); } -}; - -struct MakeMinMax { - glm::ivec2 operator()(glm::ivec3 tri) { - return glm::ivec2(glm::min(tri[0], glm::min(tri[1], tri[2])), - glm::max(tri[0], glm::max(tri[1], tri[2]))); - } -}; - -struct MinMax - : public thrust::binary_function { - glm::ivec2 operator()(glm::ivec2 a, glm::ivec2 b) { - a[0] = glm::min(a[0], b[0]); - a[1] = glm::max(a[1], b[1]); - return a; - } -}; - -struct SumPair : public thrust::binary_function, - thrust::pair, - thrust::pair> { - thrust::pair operator()(thrust::pair a, - thrust::pair b) { - a.first += b.first; - a.second += b.second; - return a; + return std::make_pair(area / 2.0f, volume / 6.0f); } }; @@ -375,10 +326,18 @@ void Manifold::Impl::CalculateBBox() { auto policy = autoPolicy(NumVert()); bBox_.min = reduce( policy, vertPos_.begin(), vertPos_.end(), - glm::vec3(std::numeric_limits::infinity()), PosMin()); + glm::vec3(std::numeric_limits::infinity()), [](auto a, auto b) { + if (isnan(a.x)) return b; + if (isnan(b.x)) return a; + return glm::min(a, b); + }); bBox_.max = reduce( policy, vertPos_.begin(), vertPos_.end(), - glm::vec3(-std::numeric_limits::infinity()), PosMax()); + glm::vec3(-std::numeric_limits::infinity()), [](auto a, auto b) { + if (isnan(a.x)) return b; + if (isnan(b.x)) return a; + return glm::max(a, b); + }); } /** @@ -387,9 +346,10 @@ void Manifold::Impl::CalculateBBox() { */ bool Manifold::Impl::IsFinite() const { auto policy = autoPolicy(NumVert()); - return transform_reduce(policy, vertPos_.begin(), vertPos_.end(), - FiniteVert(), true, - thrust::logical_and()); + return transform_reduce( + policy, vertPos_.begin(), vertPos_.end(), true, + [](bool a, bool b) { return a && b; }, + [](auto v) { return glm::all(glm::isfinite(v)); }); } /** @@ -399,10 +359,18 @@ bool Manifold::Impl::IsFinite() const { bool Manifold::Impl::IsIndexInBounds(VecView triVerts) const { auto policy = autoPolicy(triVerts.size()); glm::ivec2 minmax = transform_reduce( - policy, triVerts.begin(), triVerts.end(), MakeMinMax(), + policy, triVerts.begin(), triVerts.end(), glm::ivec2(std::numeric_limits::max(), std::numeric_limits::min()), - MinMax()); + [](auto a, auto b) { + a[0] = glm::min(a[0], b[0]); + a[1] = glm::max(a[1], b[1]); + return a; + }, + [](auto tri) { + return glm::ivec2(glm::min(tri[0], glm::min(tri[1], tri[2])), + glm::max(tri[0], glm::max(tri[1], tri[2]))); + }); return minmax[0] >= 0 && minmax[1] < static_cast(NumVert()); } @@ -430,6 +398,8 @@ float Manifold::Impl::MinGap(const Manifold::Impl& other, float minDistanceSquared = transform_reduce( autoPolicy(collisions.size()), countAt(0_z), countAt(collisions.size()), + searchLength * searchLength, + [](float a, float b) { return std::min(a, b); }, [&collisions, this, &other](int i) { const int tri = collisions.Get(i, 1); const int triOther = collisions.Get(i, 0); @@ -443,8 +413,7 @@ float Manifold::Impl::MinGap(const Manifold::Impl& other, } return DistanceTriangleTriangleSquared(p, q); - }, - searchLength * searchLength, thrust::minimum()); + }); return sqrt(minDistanceSquared); }; diff --git a/src/manifold/src/sort.cpp b/src/manifold/src/sort.cpp index 2e156d178..108631551 100644 --- a/src/manifold/src/sort.cpp +++ b/src/manifold/src/sort.cpp @@ -24,7 +24,7 @@ using namespace manifold; constexpr uint32_t kNoCode = 0xFFFFFFFFu; -struct Extrema : public thrust::binary_function { +struct Extrema { void MakeForward(Halfedge& a) { if (!a.IsForward()) { int tmp = a.startVert; @@ -112,22 +112,6 @@ struct ReindexFace { } } }; - -struct Duplicate { - thrust::pair operator()(float x) { - return thrust::make_pair(x, x); - } -}; - -struct MinMax : public thrust::binary_function, - thrust::pair, - thrust::pair> { - thrust::pair operator()(thrust::pair a, - thrust::pair b) { - return thrust::make_pair(glm::min(a.first, b.first), - glm::max(a.second, b.second)); - } -}; } // namespace namespace manifold { @@ -490,11 +474,15 @@ bool MeshGL::Merge() { Box bBox; for (const int i : {0, 1, 2}) { auto iPos = StridedRange(vertPropD.begin() + i, vertPropD.end(), numProp); - auto minMax = transform_reduce>( - autoPolicy(numVert), iPos.begin(), iPos.end(), Duplicate(), - thrust::make_pair(std::numeric_limits::infinity(), - -std::numeric_limits::infinity()), - MinMax()); + auto minMax = transform_reduce>( + autoPolicy(numVert), iPos.begin(), iPos.end(), + std::make_pair(std::numeric_limits::infinity(), + -std::numeric_limits::infinity()), + [](auto a, auto b) { + return std::make_pair(glm::min(a.first, b.first), + glm::max(a.second, b.second)); + }, + [](float f) { return std::make_pair(f, f); }); bBox.min[i] = minMax.first; bBox.max[i] = minMax.second; } diff --git a/src/polygon/src/polygon.cpp b/src/polygon/src/polygon.cpp index 2b931f1e3..2bf39d930 100644 --- a/src/polygon/src/polygon.cpp +++ b/src/polygon/src/polygon.cpp @@ -191,9 +191,9 @@ bool IsConvex(const PolygonsIdx &polys, float precision) { */ std::vector TriangulateConvex(const PolygonsIdx &polys) { const size_t numTri = transform_reduce( - autoPolicy(polys.size()), polys.begin(), polys.end(), - [](const SimplePolygonIdx &poly) { return poly.size() - 2; }, 0, - thrust::plus()); + autoPolicy(polys.size()), polys.begin(), polys.end(), 0, + [](size_t a, size_t b) { return a + b; }, + [](const SimplePolygonIdx &poly) { return poly.size() - 2; }); std::vector triangles; triangles.reserve(numTri); for (const SimplePolygonIdx &poly : polys) { diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index 561ee9360..1aeb58606 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -18,18 +18,6 @@ add_library(${PROJECT_NAME} INTERFACE) message("Parallel Backend: ${MANIFOLD_PAR}") include(FetchContent) -FetchContent_Declare(Thrust - GIT_REPOSITORY https://github.com/NVIDIA/thrust.git - GIT_TAG 2.1.0 - GIT_SHALLOW TRUE - GIT_PROGRESS TRUE -) -find_package(Thrust QUIET) -if(NOT Thrust_FOUND AND NOT DEFINED thrust_SOURCE_DIR) - message("thrust not found, downloading from source") - FetchContent_Populate(Thrust) -endif() - if (TRACY_ENABLE) include(FetchContent) FetchContent_Declare(tracy @@ -50,6 +38,13 @@ if(MANIFOLD_PAR STREQUAL "TBB") target_include_directories(${PROJECT_NAME} INTERFACE $) target_link_libraries(${PROJECT_NAME} INTERFACE ${TBB_LINK_LIBRARIES}) endif() + if(APPLE) + if(oneDPL_FOUND) + target_link_libraries(${PROJECT_NAME} INTERFACE oneDPL) + else() + message(WARNING "oneDPL not found, sequential implementation is used instead") + endif() + endif() elseif(MANIFOLD_PAR STREQUAL "NONE") set(MANIFOLD_PAR "CPP") else() @@ -62,19 +57,6 @@ target_include_directories(${PROJECT_NAME} INTERFACE $) target_link_libraries(${PROJECT_NAME} INTERFACE glm::glm) -if(NOT DEFINED thrust_SOURCE_DIR) - set(thrust_SOURCE_DIR ${_THRUST_INCLUDE_DIR}) -endif() - -target_include_directories(${PROJECT_NAME} INTERFACE - $ - $ -) - -target_compile_options(${PROJECT_NAME} INTERFACE - -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_${MANIFOLD_PAR} -) - if(MANIFOLD_EXCEPTIONS) target_compile_options(${PROJECT_NAME} INTERFACE -DMANIFOLD_EXCEPTIONS=1 diff --git a/src/utilities/include/iters.h b/src/utilities/include/iters.h index 7f7408a39..332883bfe 100644 --- a/src/utilities/include/iters.h +++ b/src/utilities/include/iters.h @@ -44,9 +44,9 @@ struct TransformIterator { F f; public: - // users are not suppposed to take pointer/reference of the iterator. using pointer = void; - using reference = void; + using reference = + std::invoke_result_t::value_type>; using difference_type = typename InnerIter::difference_type; using value_type = std::invoke_result_t::value_type>; @@ -54,6 +54,13 @@ struct TransformIterator { TransformIterator(Iter iter, F f) : iter(iter), f(f) {} + TransformIterator& operator=(const TransformIterator& other) { + if (this == &other) return *this; + // don't copy function, should be the same + iter = other.iter; + return *this; + } + value_type operator*() const { return f(*iter); } value_type operator[](size_t i) const { return f(iter[i]); } @@ -71,6 +78,19 @@ struct TransformIterator { return old; } + // prefix increment + TransformIterator& operator--() { + iter -= 1; + return *this; + } + + // postfix + TransformIterator operator--(int) { + auto old = *this; + operator--(); + return old; + } + TransformIterator operator+(size_t n) const { return TransformIterator(iter + n, f); } @@ -80,20 +100,25 @@ struct TransformIterator { return *this; } - friend bool operator==(TransformIterator a, TransformIterator b) { - return a.iter == b.iter; + TransformIterator operator-(size_t n) const { + return TransformIterator(iter - n, f); } - friend bool operator!=(TransformIterator a, TransformIterator b) { - return !(a.iter == b.iter); + TransformIterator& operator-=(size_t n) { + iter -= n; + return *this; } - friend bool operator<(TransformIterator a, TransformIterator b) { - return a.iter < b.iter; + bool operator==(TransformIterator other) const { return iter == other.iter; } + + bool operator!=(TransformIterator other) const { + return !(iter == other.iter); } - friend difference_type operator-(TransformIterator a, TransformIterator b) { - return a.iter - b.iter; + bool operator<(TransformIterator other) const { return iter < other.iter; } + + difference_type operator-(TransformIterator other) const { + return iter - other.iter; } operator TransformIterator() const { @@ -131,6 +156,19 @@ struct CountingIterator { return old; } + // prefix increment + CountingIterator& operator--() { + counter -= 1; + return *this; + } + + // postfix + CountingIterator operator--(int) { + auto old = *this; + operator--(); + return old; + } + CountingIterator operator+(T n) const { return CountingIterator(counter + n); } @@ -140,6 +178,15 @@ struct CountingIterator { return *this; } + CountingIterator operator-(T n) const { + return CountingIterator(counter - n); + } + + CountingIterator& operator-=(T n) { + counter -= n; + return *this; + } + friend bool operator==(CountingIterator a, CountingIterator b) { return a.counter == b.counter; } @@ -204,6 +251,19 @@ struct StridedRange { return old; } + // prefix increment + StridedRangeIter& operator--() { + iter -= stride; + return *this; + } + + // postfix + StridedRangeIter operator--(int) { + auto old = *this; + operator--(); + return old; + } + StridedRangeIter operator+(size_t n) const { return StridedRangeIter(iter + n * stride, stride); } @@ -213,6 +273,15 @@ struct StridedRange { return *this; } + StridedRangeIter operator-(size_t n) const { + return StridedRangeIter(iter - n * stride, stride); + } + + StridedRangeIter& operator-=(size_t n) { + iter -= n * stride; + return *this; + } + friend bool operator==(StridedRangeIter a, StridedRangeIter b) { return a.iter == b.iter; } diff --git a/src/utilities/include/par.h b/src/utilities/include/par.h index fb29b48d0..99fa82772 100644 --- a/src/utilities/include/par.h +++ b/src/utilities/include/par.h @@ -13,33 +13,23 @@ // limitations under the License. #pragma once -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include #if MANIFOLD_PAR == 'T' -#include - -#if MANIFOLD_PAR == 'T' && TBB_INTERFACE_VERSION >= 10000 && \ - __has_include() +#if __has_include() #include +#define HAS_PAR_UNSEQ +#elif __has_include() +#include +#include +#include +#include +#define HAS_PAR_UNSEQ #endif - -#define MANIFOLD_PAR_NS tbb -#else -#define MANIFOLD_PAR_NS cpp #endif +#include +#include #include "iters.h" +#include "public.h" namespace manifold { enum class ExecutionPolicy { @@ -59,34 +49,7 @@ inline constexpr ExecutionPolicy autoPolicy(size_t size) { return ExecutionPolicy::Par; } -#define THRUST_DYNAMIC_BACKEND_VOID(NAME) \ - template \ - void NAME(ExecutionPolicy policy, Args... args) { \ - switch (policy) { \ - case ExecutionPolicy::Par: \ - thrust::NAME(thrust::MANIFOLD_PAR_NS::par, args...); \ - break; \ - case ExecutionPolicy::Seq: \ - thrust::NAME(thrust::cpp::par, args...); \ - break; \ - } \ - } - -#define THRUST_DYNAMIC_BACKEND(NAME, RET) \ - template \ - Ret NAME(ExecutionPolicy policy, Args... args) { \ - switch (policy) { \ - case ExecutionPolicy::Par: \ - return thrust::NAME(thrust::MANIFOLD_PAR_NS::par, args...); \ - case ExecutionPolicy::Seq: \ - break; \ - } \ - return thrust::NAME(thrust::cpp::par, args...); \ - } - -#if MANIFOLD_PAR != 'T' || \ - (TBB_INTERFACE_VERSION >= 10000 && __has_include()) -#if MANIFOLD_PAR == 'T' +#ifdef HAS_PAR_UNSEQ #define STL_DYNAMIC_BACKEND(NAME, RET) \ template \ Ret NAME(ExecutionPolicy policy, Args... args) { \ @@ -128,87 +91,68 @@ void exclusive_scan(ExecutionPolicy policy, Args... args) { // https://github.com/llvm/llvm-project/issues/59810 std::exclusive_scan(args...); } -template -OutputIterator copy_if(ExecutionPolicy policy, InputIterator1 first, - InputIterator1 last, InputIterator2 stencil, - OutputIterator result, Predicate pred) { - if (policy == ExecutionPolicy::Seq) - return thrust::copy_if(thrust::cpp::par, first, last, stencil, result, - pred); - else - // note: this is not a typo, see - // https://github.com/NVIDIA/thrust/issues/1977 - return thrust::copy_if(first, last, stencil, result, pred); -} -template -OutputIterator copy_if(ExecutionPolicy policy, InputIterator1 first, - InputIterator1 last, OutputIterator result, - Predicate pred) { -#if MANIFOLD_PAR == 'T' - if (policy == ExecutionPolicy::Seq) - return std::copy_if(first, last, result, pred); - else - return std::copy_if(std::execution::par_unseq, first, last, result, pred); -#else - return std::copy_if(first, last, result, pred); -#endif -} - -#else -#define STL_DYNAMIC_BACKEND(NAME, RET) THRUST_DYNAMIC_BACKEND(NAME, RET) -#define STL_DYNAMIC_BACKEND_VOID(NAME) THRUST_DYNAMIC_BACKEND_VOID(NAME) -THRUST_DYNAMIC_BACKEND_VOID(exclusive_scan) -THRUST_DYNAMIC_BACKEND(copy_if, void) -#endif - -THRUST_DYNAMIC_BACKEND_VOID(gather) -THRUST_DYNAMIC_BACKEND_VOID(scatter) -THRUST_DYNAMIC_BACKEND_VOID(for_each) -THRUST_DYNAMIC_BACKEND_VOID(for_each_n) -THRUST_DYNAMIC_BACKEND_VOID(sequence) -STL_DYNAMIC_BACKEND_VOID(transform) -STL_DYNAMIC_BACKEND_VOID(uninitialized_fill) -STL_DYNAMIC_BACKEND_VOID(uninitialized_copy) -STL_DYNAMIC_BACKEND_VOID(stable_sort) -STL_DYNAMIC_BACKEND_VOID(fill) -STL_DYNAMIC_BACKEND_VOID(inclusive_scan) +template +void scatter(ExecutionPolicy policy, InputIterator1 first, InputIterator1 last, + InputIterator2 mapFirst, OutputIterator outputFirst) { + for_each(policy, countAt(0_z), + countAt(static_cast(std::distance(first, last))), + [first, mapFirst, outputFirst](size_t i) { + outputFirst[mapFirst[i]] = first[i]; + }); +} -// there are some issues with thrust copy -template -OutputIterator copy(ExecutionPolicy policy, InputIterator first, - InputIterator last, OutputIterator result) { -#if MANIFOLD_PAR == 'T' && \ - (TBB_INTERFACE_VERSION >= 10000 && __has_include()) - if (policy == ExecutionPolicy::Par) - return std::copy(std::execution::par_unseq, first, last, result); -#endif - return std::copy(first, last, result); +template +void gather(ExecutionPolicy policy, InputIterator mapFirst, + InputIterator mapLast, RandomAccessIterator inputFirst, + OutputIterator outputFirst) { + for_each(policy, countAt(0_z), + countAt(static_cast(std::distance(mapFirst, mapLast))), + [mapFirst, inputFirst, outputFirst](size_t i) { + outputFirst[i] = inputFirst[mapFirst[i]]; + }); } -template -OutputIterator copy_n(ExecutionPolicy policy, InputIterator first, size_t n, - OutputIterator result) { -#if MANIFOLD_PAR == 'T' && \ - (TBB_INTERFACE_VERSION >= 10000 && __has_include()) - if (policy == ExecutionPolicy::Par) - return std::copy_n(std::execution::par_unseq, first, n, result); -#endif - return std::copy_n(first, n, result); +template +void sequence(ExecutionPolicy policy, Iterator first, Iterator last) { + for_each(policy, countAt(0_z), + countAt(static_cast(std::distance(first, last))), + [first](size_t i) { first[i] = i; }); } // void implies that the user have to specify the return type in the template // argument, as we are unable to deduce it -THRUST_DYNAMIC_BACKEND(transform_reduce, void) STL_DYNAMIC_BACKEND(remove, void) STL_DYNAMIC_BACKEND(find, void) STL_DYNAMIC_BACKEND(find_if, void) STL_DYNAMIC_BACKEND(all_of, bool) STL_DYNAMIC_BACKEND(is_sorted, bool) -STL_DYNAMIC_BACKEND(reduce, void) +// STL_DYNAMIC_BACKEND(reduce, void) +template +Ret reduce(ExecutionPolicy policy, Args... args) { + return std::reduce(args...); +} STL_DYNAMIC_BACKEND(count_if, int) STL_DYNAMIC_BACKEND(remove_if, void) +STL_DYNAMIC_BACKEND(copy_if, void) +STL_DYNAMIC_BACKEND(unique, void) +// STL_DYNAMIC_BACKEND(transform_reduce, void) +template +Ret transform_reduce(ExecutionPolicy policy, Args... args) { + return std::transform_reduce(args...); +} + +STL_DYNAMIC_BACKEND_VOID(for_each) +STL_DYNAMIC_BACKEND_VOID(for_each_n) +STL_DYNAMIC_BACKEND_VOID(transform) +STL_DYNAMIC_BACKEND_VOID(uninitialized_fill) +STL_DYNAMIC_BACKEND_VOID(uninitialized_copy) +STL_DYNAMIC_BACKEND_VOID(stable_sort) +STL_DYNAMIC_BACKEND_VOID(fill) +STL_DYNAMIC_BACKEND_VOID(inclusive_scan) +STL_DYNAMIC_BACKEND_VOID(copy) +STL_DYNAMIC_BACKEND_VOID(copy_n) } // namespace manifold diff --git a/src/utilities/include/sparse.h b/src/utilities/include/sparse.h index 9bd1fc8e6..9615f10f2 100644 --- a/src/utilities/include/sparse.h +++ b/src/utilities/include/sparse.h @@ -128,7 +128,9 @@ class SparseIndices { void Unique() { Sort(); VecView view = AsVec64(); - size_t newSize = std::unique(view.begin(), view.end()) - view.begin(); + size_t newSize = unique(autoPolicy(view.size()), + view.begin(), view.end()) - + view.begin(); Resize(newSize); } diff --git a/src/utilities/include/utils.h b/src/utilities/include/utils.h index e6b3c41b7..ba81c8407 100644 --- a/src/utilities/include/utils.h +++ b/src/utilities/include/utils.h @@ -184,5 +184,16 @@ struct UnionFind { return toLabel.size() + lonelyNodes; } }; + +template +struct Identity { + T operator()(T v) const { return v; } +}; + +template +struct Negate { + T operator()(T v) const { return -v; } +}; + /** @} */ } // namespace manifold