OpenMP Parallelization for CascadedUnion [WIP]

CMakeLists.txt

@@ -174,6 +174,7 @@ target_compile_options(geos_cxx_flags INTERFACE
   $<$<AND:$<CXX_COMPILER_ID:GNU>,$<EQUAL:4,${CMAKE_SIZEOF_VOID_P}>>:-ffloat-store>
 )
 
+
 #-----------------------------------------------------------------------------
 # Target geos_cxx_flags: common compilation flags
 #-----------------------------------------------------------------------------
@@ -260,6 +261,31 @@ target_link_libraries(geos PUBLIC geos_cxx_flags PRIVATE $<BUILD_INTERFACE:ryu>)
 add_subdirectory(include)
 add_subdirectory(src)
 
+#-----------------------------------------------------------------------------
+# OpenMP
+#-----------------------------------------------------------------------------
+
+option(GEOS_ENABLE_OPENMP "Enable OpenMP" ON)
+if(GEOS_ENABLE_OPENMP)
+  find_package(OpenMP)
+  if(OpenMP_CXX_FOUND)
+    message(STATUS "GEOS: OpenMP_C_INCLUDE_DIRS ${OpenMP_C_INCLUDE_DIRS}")
+    message(STATUS "GEOS: OpenMP_CXX_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS}")
+    message(STATUS "GEOS: OpenMP_CXX_LIBRARY ${OpenMP_CXX_LIBRARY}")
+    message(STATUS "GEOS: OpenMP_CXX_LIBRARIES ${OpenMP_CXX_LIBRARIES}")
+    message(STATUS "GEOS: OpenMP_CXX_LIB_NAMES ${OpenMP_CXX_LIB_NAMES}")
+    message(STATUS "GEOS: OpenMP_CXX_FLAGS ${OpenMP_CXX_FLAGS}")
+    target_compile_options(geos_cxx_flags INTERFACE "${OpenMP_CXX_FLAGS}")
+    target_compile_definitions(geos_cxx_flags INTERFACE GEOS_OPENMP)
+    target_include_directories(geos INTERFACE ${OpenMP_CXX_INCLUDE_DIRS})
+    target_link_libraries(geos PRIVATE ${OpenMP_CXX_LIBRARIES})
+  else()
+    message(STATUS "GEOS: OpenMP not found")
+  endif()
+endif()
+
+#-----------------------------------------------------------------------------
+
 if(BUILD_SHARED_LIBS)
   target_compile_definitions(geos
     PRIVATE $<IF:$<CXX_COMPILER_ID:MSVC>,GEOS_DLL_EXPORT,DLL_EXPORT>)

include/geos/operation/union/CascadedPolygonUnion.h

@@ -99,10 +99,21 @@ class GEOS_DLL ClassicUnionStrategy : public UnionStrategy {
  * between the input geometries. However, this case is likely rare in practice.
  */
 class GEOS_DLL CascadedPolygonUnion {
+
 private:
+
     std::vector<geom::Polygon*>* inputPolys;
     geom::GeometryFactory const* geomFactory;
 
+    /* Default when not called from SR or NG modules */
+    UnionStrategy* unionFunction;
+    ClassicUnionStrategy defaultUnionFunction;
+
+    /* private */
+    std::unique_ptr<geom::Geometry>
+
+    pairUnion(std::vector<const geom::Geometry*>& polys) const;
+
     /**
      * The effectiveness of the index is somewhat sensitive
      * to the node capacity.
@@ -126,10 +137,50 @@ class GEOS_DLL CascadedPolygonUnion {
      * @param g the geometry to filter
      * @return a Polygonal geometry
      */
-    static std::unique_ptr<geom::Geometry> restrictToPolygons(std::unique_ptr<geom::Geometry> g);
+    static std::unique_ptr<geom::Geometry> restrictToPolygons(
+        std::unique_ptr<geom::Geometry> g);
+
+    /**
+     * Unions a section of a list using a recursive binary union on each half
+     * of the section.
+     *
+     * @param geoms the list of geometries containing the section to union
+     * @param start the start index of the section
+     * @param end the index after the end of the section
+     * @return the union of the list section
+     */
+    std::unique_ptr<geom::Geometry> binaryUnion(
+        const std::vector<const geom::Geometry*> & geoms,
+        std::size_t start, std::size_t end) const;
+
+    /**
+     * Computes the union of two geometries,
+     * either of both of which may be null.
+     *
+     * @param g0 a geom::Geometry
+     * @param g1 a geom::Geometry
+     * @return the union of the input(s) or null if both inputs are null
+     */
+    std::unique_ptr<geom::Geometry> unionSafe(
+        const geom::Geometry* g0, const geom::Geometry* g1) const;
+
+    std::unique_ptr<geom::Geometry> unionSafe(
+        std::unique_ptr<geom::Geometry> &&, std::unique_ptr<geom::Geometry> &&) const;
+
+    /**
+     * Encapsulates the actual unioning of two polygonal geometries.
+     *
+     * @param g0 a geom::Geometry
+     * @param g1 a geom::Geometry
+     * @return The union of the inputs
+     */
+    std::unique_ptr<geom::Geometry> unionActual(
+        const geom::Geometry* g0, const geom::Geometry* g1) const;
+
+    std::unique_ptr<geom::Geometry> unionActual(
+        std::unique_ptr<geom::Geometry> &&, std::unique_ptr<geom::Geometry> &&) const;
 
 public:
-    CascadedPolygonUnion();
 
     /** \brief
      * Computes the union of a collection of polygonal [Geometrys](@ref geom::Geometry).
@@ -172,68 +223,29 @@ class GEOS_DLL CascadedPolygonUnion {
      * Creates a new instance to union the given collection of
      * [Geometrys](@ref geom::Geometry).
      *
-     * @param polys a collection of polygonal [Geometrys](@ref geom::Geometry).
-     *              Ownership of elements *and* vector are left to caller.
+     * @param polys A collection of polygonal [Geometrys](@ref geom::Geometry).
+     *        Ownership of elements *and* vector are left to caller.
+     * @param unionFun Function to call for pairwise unions. Might implement
+     *        an alternate strategy (snap rounding vs snapping vs full precision)
      */
-    CascadedPolygonUnion(std::vector<geom::Polygon*>* polys)
-        : inputPolys(polys)
-        , geomFactory(nullptr)
-        , unionFunction(&defaultUnionFunction)
-    {}
-
     CascadedPolygonUnion(std::vector<geom::Polygon*>* polys, UnionStrategy* unionFun)
         : inputPolys(polys)
         , geomFactory(nullptr)
         , unionFunction(unionFun)
     {}
 
+    CascadedPolygonUnion(std::vector<geom::Polygon*>* polys)
+        : CascadedPolygonUnion(polys, &defaultUnionFunction)
+    {}
+
     /** \brief
      * Computes the union of the input geometries.
      *
-     * @return the union of the input geometries
-     * @return `null` if no input geometries were provided
-     */
-    std::unique_ptr<geom::Geometry> Union();
-
-private:
-
-    UnionStrategy* unionFunction;
-    ClassicUnionStrategy defaultUnionFunction;
-
-    /**
-     * Unions a section of a list using a recursive binary union on each half
-     * of the section.
-     *
-     * @param geoms the list of geometries containing the section to union
-     * @param start the start index of the section
-     * @param end the index after the end of the section
-     * @return the union of the list section
-     */
-    std::unique_ptr<geom::Geometry> binaryUnion(const std::vector<const geom::Geometry*> & geoms, std::size_t start, std::size_t end);
-
-    /**
-     * Computes the union of two geometries,
-     * either of both of which may be null.
-     *
-     * @param g0 a Geometry
-     * @param g1 a Geometry
-     * @return the union of the input(s)
-     * @return null if both inputs are null
-     */
-    std::unique_ptr<geom::Geometry> unionSafe(const geom::Geometry* g0, const geom::Geometry* g1) const;
-
-    std::unique_ptr<geom::Geometry> unionSafe(std::unique_ptr<geom::Geometry> &&, std::unique_ptr<geom::Geometry> &&);
-
-    /**
-     * Encapsulates the actual unioning of two polygonal geometries.
+     * @return the union of the input geometries or `nullptr` if no input geometries were provided
      *
-     * @param g0
-     * @param g1
-     * @return
      */
-    std::unique_ptr<geom::Geometry> unionActual(const geom::Geometry* g0, const geom::Geometry* g1) const;
+    std::unique_ptr<geom::Geometry> Union();
 
-    std::unique_ptr<geom::Geometry> unionActual(std::unique_ptr<geom::Geometry> &&, std::unique_ptr<geom::Geometry> &&) const;
 };
 
 

src/operation/union/CascadedPolygonUnion.cpp

@@ -42,21 +42,23 @@ namespace operation { // geos.operation
 namespace geounion {  // geos.operation.geounion
 
 
-// ////////////////////////////////////////////////////////////////////////////
+/* public static */
 std::unique_ptr<geom::Geometry>
 CascadedPolygonUnion::Union(std::vector<geom::Polygon*>* polys)
 {
     CascadedPolygonUnion op(polys);
     return op.Union();
 }
 
+/* public static */
 std::unique_ptr<geom::Geometry>
 CascadedPolygonUnion::Union(std::vector<geom::Polygon*>* polys, UnionStrategy* unionFun)
 {
     CascadedPolygonUnion op(polys, unionFun);
     return op.Union();
 }
 
+/* public static */
 std::unique_ptr<geom::Geometry>
 CascadedPolygonUnion::Union(const geom::MultiPolygon* multipoly)
 {
@@ -70,6 +72,7 @@ CascadedPolygonUnion::Union(const geom::MultiPolygon* multipoly)
     return op.Union();
 }
 
+/* public */
 std::unique_ptr<geom::Geometry>
 CascadedPolygonUnion::Union()
 {
@@ -94,13 +97,61 @@ CascadedPolygonUnion::Union()
     // TODO avoid creating this vector and run binaryUnion off the iterators directly
     std::vector<const geom::Geometry*> geoms(index.items().begin(), index.items().end());
 
-    return binaryUnion(geoms, 0, geoms.size());
+// #ifdef GEOS_OPENMP
+    return pairUnion(geoms);
+// #else
+//     return binaryUnion(geoms, 0, geoms.size());
+// #endif
 }
 
+#ifdef GEOS_OPENMP
+#include <omp.h>
+#endif
 
+/* private */
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::binaryUnion(const std::vector<const geom::Geometry*> & geoms,
-                                  std::size_t start, std::size_t end)
+CascadedPolygonUnion::pairUnion(
+    std::vector<const geom::Geometry*>& inply) const
+{
+    /* Windows OpenMP requires loop iterators to be int */
+    /* so we do a lot of dumb casting up here to allow ints */
+    /* in the loops further down */
+    int sz = (int)(inply.size());
+    int nsz = (sz / 2) + (sz % 2);
+    std::vector<std::unique_ptr<geom::Geometry>> outply((std::size_t)nsz);
+
+#pragma omp parallel for
+    for (int i = 0; i < sz; i += 2) {
+        /* Work backwards to preserve an old regression result */
+        const geom::Geometry* g0 = inply[(std::size_t)(sz - i - 1)];
+        const geom::Geometry* g1 = i+1 < sz ? inply[(std::size_t)(sz - i - 2)] : nullptr;
+        outply[(std::size_t)(i/2)] = unionSafe(g0, g1);
+    }
+
+    int step = 1;
+    while(step < nsz) {
+
+#pragma omp parallel for
+        /* Place resultant output into position 1 of the inputs */
+        /* This keeps threads from stomping on each other as they */
+        /* work through different parts of the processing */
+        for (int i0 = 0; i0 < nsz; i0 += 2*step) {
+            int i1 = i0 + step;
+            const geom::Geometry* g0 = outply[(std::size_t)(i0)].release();
+            const geom::Geometry* g1 = i1 < nsz ? outply[(std::size_t)(i1)].release() : nullptr;
+            outply[(std::size_t)(i0)] = unionSafe(g0, g1);
+        }
+        step *= 2;
+    }
+    return std::move(outply[0]);
+}
+
+
+/* private */
+std::unique_ptr<geom::Geometry>
+CascadedPolygonUnion::binaryUnion(
+    const std::vector<const geom::Geometry*> & geoms,
+    std::size_t start, std::size_t end) const
 {
     if(end - start <= 1) {
         return unionSafe(geoms[start], nullptr);
@@ -117,8 +168,10 @@ CascadedPolygonUnion::binaryUnion(const std::vector<const geom::Geometry*> & geo
     }
 }
 
+/* private */
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::unionSafe(const geom::Geometry* g0, const geom::Geometry* g1) const
+CascadedPolygonUnion::unionSafe(
+    const geom::Geometry* g0, const geom::Geometry* g1) const
 {
     if(g0 == nullptr && g1 == nullptr) {
         return nullptr;
@@ -134,8 +187,11 @@ CascadedPolygonUnion::unionSafe(const geom::Geometry* g0, const geom::Geometry*
     return unionActual(g0, g1);
 }
 
+/* private */
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::unionSafe(std::unique_ptr<geom::Geometry> && g0, std::unique_ptr<geom::Geometry> && g1)
+CascadedPolygonUnion::unionSafe(
+    std::unique_ptr<geom::Geometry> && g0,
+    std::unique_ptr<geom::Geometry> && g1) const
 {
     if(g0 == nullptr && g1 == nullptr) {
         return nullptr;
@@ -151,22 +207,29 @@ CascadedPolygonUnion::unionSafe(std::unique_ptr<geom::Geometry> && g0, std::uniq
     return unionActual(std::move(g0), std::move(g1));
 }
 
+/* private */
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::unionActual(const geom::Geometry* g0, const geom::Geometry* g1) const
+CascadedPolygonUnion::unionActual(
+    const geom::Geometry* g0,
+    const geom::Geometry* g1) const
 {
     std::unique_ptr<geom::Geometry> ug;
     ug = unionFunction->Union(g0, g1);
     return restrictToPolygons(std::move(ug));
 }
 
+/* private */
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::unionActual(std::unique_ptr<geom::Geometry> && g0, std::unique_ptr<geom::Geometry> && g1) const
+CascadedPolygonUnion::unionActual(
+    std::unique_ptr<geom::Geometry> && g0,
+    std::unique_ptr<geom::Geometry> && g1) const
 {
     std::unique_ptr<geom::Geometry> ug;
     ug = unionFunction->Union(std::move(g0), std::move(g1));
     return restrictToPolygons(std::move(ug));
 }
 
+/* private static */
 std::unique_ptr<geom::Geometry>
 CascadedPolygonUnion::restrictToPolygons(std::unique_ptr<geom::Geometry> g)
 {