From 92a18b4b144ffadd9d884065a156bc948c1a310b Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Fri, 16 Jul 2021 13:49:44 -0700 Subject: [PATCH 1/4] Parallel processing of cascaded union using OpenMP --- CMakeLists.txt | 23 ++++ .../operation/union/CascadedPolygonUnion.h | 114 ++++++++++-------- src/operation/union/CascadedPolygonUnion.cpp | 73 +++++++++-- 3 files changed, 151 insertions(+), 59 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e99aa0a171..e6b5ae49fe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -174,6 +174,7 @@ target_compile_options(geos_cxx_flags INTERFACE $<$,$>:-ffloat-store> ) + #----------------------------------------------------------------------------- # Target geos_cxx_flags: common compilation flags #----------------------------------------------------------------------------- @@ -260,6 +261,28 @@ target_link_libraries(geos PUBLIC geos_cxx_flags PRIVATE $) 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_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_include_directories(geos INTERFACE ${OpenMP_CXX_INCLUDE_DIRS}) + # target_link_options(geos PUBLIC "LINKER:-L${OpenMP_CXX_LIBRARY},-l${OpenMP_CXX_LIB_NAMES}") + target_compile_definitions(geos_cxx_flags INTERFACE GEOS_OPENMP) + endif() +endif() + +#----------------------------------------------------------------------------- + if(BUILD_SHARED_LIBS) target_compile_definitions(geos PRIVATE $,GEOS_DLL_EXPORT,DLL_EXPORT>) diff --git a/include/geos/operation/union/CascadedPolygonUnion.h b/include/geos/operation/union/CascadedPolygonUnion.h index 0bcf39426f..3fe17219fd 100644 --- a/include/geos/operation/union/CascadedPolygonUnion.h +++ b/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* inputPolys; geom::GeometryFactory const* geomFactory; + /* Default when not called from SR or NG modules */ + UnionStrategy* unionFunction; + ClassicUnionStrategy defaultUnionFunction; + + /* private */ + std::unique_ptr + + pairUnion(std::vector& 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 restrictToPolygons(std::unique_ptr g); + static std::unique_ptr restrictToPolygons( + std::unique_ptr 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 binaryUnion( + const std::vector & 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 unionSafe( + const geom::Geometry* g0, const geom::Geometry* g1) const; + + std::unique_ptr unionSafe( + std::unique_ptr &&, std::unique_ptr &&) 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 unionActual( + const geom::Geometry* g0, const geom::Geometry* g1) const; + + std::unique_ptr unionActual( + std::unique_ptr &&, std::unique_ptr &&) 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* polys) - : inputPolys(polys) - , geomFactory(nullptr) - , unionFunction(&defaultUnionFunction) - {} - CascadedPolygonUnion(std::vector* polys, UnionStrategy* unionFun) : inputPolys(polys) , geomFactory(nullptr) , unionFunction(unionFun) {} + CascadedPolygonUnion(std::vector* 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 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 binaryUnion(const std::vector & 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 unionSafe(const geom::Geometry* g0, const geom::Geometry* g1) const; - - std::unique_ptr unionSafe(std::unique_ptr &&, std::unique_ptr &&); - - /** - * 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 unionActual(const geom::Geometry* g0, const geom::Geometry* g1) const; + std::unique_ptr Union(); - std::unique_ptr unionActual(std::unique_ptr &&, std::unique_ptr &&) const; }; diff --git a/src/operation/union/CascadedPolygonUnion.cpp b/src/operation/union/CascadedPolygonUnion.cpp index c30c6f7eb5..db4ceb2903 100644 --- a/src/operation/union/CascadedPolygonUnion.cpp +++ b/src/operation/union/CascadedPolygonUnion.cpp @@ -42,7 +42,7 @@ namespace operation { // geos.operation namespace geounion { // geos.operation.geounion -// //////////////////////////////////////////////////////////////////////////// +/* public static */ std::unique_ptr CascadedPolygonUnion::Union(std::vector* polys) { @@ -50,6 +50,7 @@ CascadedPolygonUnion::Union(std::vector* polys) return op.Union(); } +/* public static */ std::unique_ptr CascadedPolygonUnion::Union(std::vector* polys, UnionStrategy* unionFun) { @@ -57,6 +58,7 @@ CascadedPolygonUnion::Union(std::vector* polys, UnionStrategy* u return op.Union(); } +/* public static */ std::unique_ptr CascadedPolygonUnion::Union(const geom::MultiPolygon* multipoly) { @@ -70,6 +72,7 @@ CascadedPolygonUnion::Union(const geom::MultiPolygon* multipoly) return op.Union(); } +/* public */ std::unique_ptr CascadedPolygonUnion::Union() { @@ -94,13 +97,55 @@ CascadedPolygonUnion::Union() // TODO avoid creating this vector and run binaryUnion off the iterators directly std::vector geoms(index.items().begin(), index.items().end()); - return binaryUnion(geoms, 0, geoms.size()); + // return binaryUnion(geoms, 0, geoms.size()); + return pairUnion(geoms); } +#ifdef GEOS_OPENMP +#include +#endif +/* private */ std::unique_ptr -CascadedPolygonUnion::binaryUnion(const std::vector & geoms, - std::size_t start, std::size_t end) +CascadedPolygonUnion::pairUnion( + std::vector& inply) const +{ + std::size_t sz = inply.size(); + std::size_t nsz = (sz / 2) + (sz % 2); + std::vector> outply(nsz); + +#pragma omp parallel for + for (std::size_t i = 0; i < sz; i += 2) { + /* Work backwards to preserve an old regression result */ + const geom::Geometry* g0 = inply[sz - i - 1]; + const geom::Geometry* g1 = i+1 < sz ? inply[sz - i - 2] : nullptr; + outply[i/2] = unionSafe(g0, g1); + } + + std::size_t 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 (std::size_t i0 = 0; i0 < sz; i0 += 2*step) { + std::size_t i1 = i0 + step; + const geom::Geometry* g0 = outply[i0].release(); + const geom::Geometry* g1 = i1 < sz ? outply[i1].release() : nullptr; + outply[i0] = unionSafe(g0, g1); + } + step *= 2; + } + return std::move(outply[0]); +} + + +/* private */ +std::unique_ptr +CascadedPolygonUnion::binaryUnion( + const std::vector & geoms, + std::size_t start, std::size_t end) const { if(end - start <= 1) { return unionSafe(geoms[start], nullptr); @@ -117,8 +162,10 @@ CascadedPolygonUnion::binaryUnion(const std::vector & geo } } +/* private */ std::unique_ptr -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 +181,11 @@ CascadedPolygonUnion::unionSafe(const geom::Geometry* g0, const geom::Geometry* return unionActual(g0, g1); } +/* private */ std::unique_ptr -CascadedPolygonUnion::unionSafe(std::unique_ptr && g0, std::unique_ptr && g1) +CascadedPolygonUnion::unionSafe( + std::unique_ptr && g0, + std::unique_ptr && g1) const { if(g0 == nullptr && g1 == nullptr) { return nullptr; @@ -151,22 +201,29 @@ CascadedPolygonUnion::unionSafe(std::unique_ptr && g0, std::uniq return unionActual(std::move(g0), std::move(g1)); } +/* private */ std::unique_ptr -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 ug; ug = unionFunction->Union(g0, g1); return restrictToPolygons(std::move(ug)); } +/* private */ std::unique_ptr -CascadedPolygonUnion::unionActual(std::unique_ptr && g0, std::unique_ptr && g1) const +CascadedPolygonUnion::unionActual( + std::unique_ptr && g0, + std::unique_ptr && g1) const { std::unique_ptr ug; ug = unionFunction->Union(std::move(g0), std::move(g1)); return restrictToPolygons(std::move(ug)); } +/* private static */ std::unique_ptr CascadedPolygonUnion::restrictToPolygons(std::unique_ptr g) { From 7b2f6430b81962b90a9d774db26df1ce1f28bb82 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Fri, 16 Jul 2021 14:44:19 -0700 Subject: [PATCH 2/4] Ug, set the OMP iterator to signed for MSVC then cast it back to size_t everywhere for clang --- CMakeLists.txt | 7 +++-- src/operation/union/CascadedPolygonUnion.cpp | 31 +++++++++++--------- 2 files changed, 21 insertions(+), 17 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e6b5ae49fe..2e949a5fa3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -266,18 +266,19 @@ add_subdirectory(src) #----------------------------------------------------------------------------- option(GEOS_ENABLE_OPENMP "Enable OpenMP" ON) -if (GEOS_ENABLE_OPENMP) +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_include_directories(geos INTERFACE ${OpenMP_CXX_INCLUDE_DIRS}) - # target_link_options(geos PUBLIC "LINKER:-L${OpenMP_CXX_LIBRARY},-l${OpenMP_CXX_LIB_NAMES}") 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}) endif() endif() diff --git a/src/operation/union/CascadedPolygonUnion.cpp b/src/operation/union/CascadedPolygonUnion.cpp index db4ceb2903..506c7d05ab 100644 --- a/src/operation/union/CascadedPolygonUnion.cpp +++ b/src/operation/union/CascadedPolygonUnion.cpp @@ -97,8 +97,11 @@ CascadedPolygonUnion::Union() // TODO avoid creating this vector and run binaryUnion off the iterators directly std::vector 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 @@ -110,30 +113,30 @@ std::unique_ptr CascadedPolygonUnion::pairUnion( std::vector& inply) const { - std::size_t sz = inply.size(); - std::size_t nsz = (sz / 2) + (sz % 2); - std::vector> outply(nsz); + int sz = static_cast(inply.size()); + int nsz = (sz / 2) + (sz % 2); + std::vector> outply(static_cast(nsz)); #pragma omp parallel for - for (std::size_t i = 0; i < sz; i += 2) { + for (int i = 0; i < sz; i += 2) { /* Work backwards to preserve an old regression result */ - const geom::Geometry* g0 = inply[sz - i - 1]; - const geom::Geometry* g1 = i+1 < sz ? inply[sz - i - 2] : nullptr; - outply[i/2] = unionSafe(g0, g1); + const geom::Geometry* g0 = inply[static_cast(sz - i - 1)]; + const geom::Geometry* g1 = i+1 < sz ? inply[static_cast(sz - i - 2)] : nullptr; + outply[static_cast(i/2)] = unionSafe(g0, g1); } - std::size_t step = 1; + 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 (std::size_t i0 = 0; i0 < sz; i0 += 2*step) { - std::size_t i1 = i0 + step; - const geom::Geometry* g0 = outply[i0].release(); - const geom::Geometry* g1 = i1 < sz ? outply[i1].release() : nullptr; - outply[i0] = unionSafe(g0, g1); + for (int i0 = 0; i0 < sz; i0 += 2*step) { + int i1 = i0 + step; + const geom::Geometry* g0 = outply[static_cast(i0)].release(); + const geom::Geometry* g1 = i1 < sz ? outply[static_cast(i1)].release() : nullptr; + outply[static_cast(i0)] = unionSafe(g0, g1); } step *= 2; } From 2e84518c69bcf8f1710946918d13982f60608b90 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Fri, 16 Jul 2021 16:51:05 -0700 Subject: [PATCH 3/4] Fix off-by-something error --- CMakeLists.txt | 2 ++ src/operation/union/CascadedPolygonUnion.cpp | 26 ++++++++++---------- 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e949a5fa3..ec61b1d5f7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -279,6 +279,8 @@ if(GEOS_ENABLE_OPENMP) 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() diff --git a/src/operation/union/CascadedPolygonUnion.cpp b/src/operation/union/CascadedPolygonUnion.cpp index 506c7d05ab..1c566f6a60 100644 --- a/src/operation/union/CascadedPolygonUnion.cpp +++ b/src/operation/union/CascadedPolygonUnion.cpp @@ -97,11 +97,11 @@ CascadedPolygonUnion::Union() // TODO avoid creating this vector and run binaryUnion off the iterators directly std::vector geoms(index.items().begin(), index.items().end()); -#ifdef GEOS_OPENMP +// #ifdef GEOS_OPENMP return pairUnion(geoms); -#else - return binaryUnion(geoms, 0, geoms.size()); -#endif +// #else +// return binaryUnion(geoms, 0, geoms.size()); +// #endif } #ifdef GEOS_OPENMP @@ -113,16 +113,16 @@ std::unique_ptr CascadedPolygonUnion::pairUnion( std::vector& inply) const { - int sz = static_cast(inply.size()); + int sz = (int)(inply.size()); int nsz = (sz / 2) + (sz % 2); - std::vector> outply(static_cast(nsz)); + std::vector> 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[static_cast(sz - i - 1)]; - const geom::Geometry* g1 = i+1 < sz ? inply[static_cast(sz - i - 2)] : nullptr; - outply[static_cast(i/2)] = unionSafe(g0, g1); + 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; @@ -132,11 +132,11 @@ CascadedPolygonUnion::pairUnion( /* 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 < sz; i0 += 2*step) { + for (int i0 = 0; i0 < nsz; i0 += 2*step) { int i1 = i0 + step; - const geom::Geometry* g0 = outply[static_cast(i0)].release(); - const geom::Geometry* g1 = i1 < sz ? outply[static_cast(i1)].release() : nullptr; - outply[static_cast(i0)] = unionSafe(g0, g1); + 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; } From 70bf61f27d1eba0a8dc9f2370e19b57b8a6a5360 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Tue, 27 Jul 2021 14:25:27 -0700 Subject: [PATCH 4/4] Add comment note about Windows OpenMP and int for loops --- src/operation/union/CascadedPolygonUnion.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/operation/union/CascadedPolygonUnion.cpp b/src/operation/union/CascadedPolygonUnion.cpp index 1c566f6a60..a35dfbbc9c 100644 --- a/src/operation/union/CascadedPolygonUnion.cpp +++ b/src/operation/union/CascadedPolygonUnion.cpp @@ -113,6 +113,9 @@ std::unique_ptr CascadedPolygonUnion::pairUnion( std::vector& 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> outply((std::size_t)nsz);