From 3fee2b3709496fc2f0d87c7bb9e4cade31f91182 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 29 Jun 2026 06:09:28 +0100 Subject: [PATCH 1/4] Filtered_kernel: Optimisation of side_of_orienyed_sphere_3 --- .../benchmark/Filtered_kernel/CMakeLists.txt | 25 ++---- .../side_of_oriented_sphere_3.cpp | 77 +++++++++++++++++++ .../Side_of_oriented_sphere_3.h | 8 +- 3 files changed, 89 insertions(+), 21 deletions(-) create mode 100644 Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp diff --git a/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt b/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt index 290027a019c9..f5b0876f25f5 100644 --- a/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt +++ b/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt @@ -1,24 +1,11 @@ -# Created by the script cgal_create_cmake_script -# This is the CMake script for compiling a CGAL application. +# Created by the script cgal_create_CMakeLists +# This is the CMake script for compiling a set of CGAL applications. cmake_minimum_required(VERSION 3.12...3.31) -project(Filtered_kernel_test) -add_executable(bench_simple_comparisons bench_simple_comparisons.cpp) +project(Filtered_kernel_3_Benchmarks) -find_package(CGAL REQUIRED COMPONENTS Core) +# CGAL and its components +find_package(CGAL REQUIRED) -add_executable(bench_orientation_3 "orientation_3.cpp") -target_link_libraries(bench_orientation_3 CGAL::CGAL_Core) - -add_executable(bench_comparisons "orientation_3.cpp") -target_link_libraries(bench_comparisons CGAL::CGAL_Core) -set_property( - TARGET bench_comparisons - APPEND - PROPERTY COMPILE_DEFINITIONS ONLY_TEST_COMPARISONS) - -get_property( - DEF - TARGET bench_comparisons - PROPERTY COMPILE_DEFINITIONS) +create_single_source_cgal_program("side_of_oriented_sphere_3.cpp") diff --git a/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp new file mode 100644 index 000000000000..50d83ad74bdc --- /dev/null +++ b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp @@ -0,0 +1,77 @@ +// #define CGAL_SMALL_UNORDERED_MAP_STATS +// #define CGAL_PROFILE +//#define CGAL_USE_SSE2_FABS +//#define CGAL_USE_SSE2_MAX +//#define CGAL_MSVC_USE_STD_FABS // use this one with precise +#define CGAL_NDEBUG 1 +#define NDEBUG 1 + + +#include + + +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +typedef K::Point_3 Point_3; +typedef CGAL::Timer Timer; + + +int main(int argc, char* argv[]) +{ + std::locale loc = std::locale() + .combine>(std::locale("en_US.UTF8")); + std::cout.imbue(loc); + + int M = 10; // outer loop counter + int Q = 10000; // number of consecutive queries + + + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/ocean_r.xyz"); + if(argc > 2) { + auto M_ = std::atoi(argv[2]); + if(M_ <= 0) { + std::cerr << "Invalid number of iterations: " << M_ << ". Using default value of " << M << std::endl; + } else { + M = M_; + } + } + if(argc > 3) { + auto Q_ = std::atoi(argv[3]); + if(Q <= 0) { + std::cerr << "Invalid number of iterations: " << Q_ << ". Using default value of " << Q << std::endl; + } else { + Q = Q_; + } + } + std::ifstream in(filename.c_str()); + std::vector points; + Point_3 p, q; + + while(in >> p ){ + points.push_back(p); + } + + std::cout << points.size() << " points read\n"; + + std::cout << "We run " << M << " times: All four consecutive points combined with the next " << Q << " points as query" << std::endl; + Timer timer; + timer.start(); + std::size_t inside = 0; + std::size_t bound = (std::min)(points.size(), std::size_t(Q)); + for(int k = 0; k < M; k++) + for(int i = 0; i < points.size()-6; ++i) + for(int j = i; j < bound ; ++j ) + if(side_of_oriented_sphere(points[i],points[i+1], points[i+2], points[i+3], points[j]) == CGAL::ON_NEGATIVE_SIDE) ++inside; + + std ::cout << inside << std::endl; + + timer.stop(); + std::cout << "Time elapsed: " << timer.time() << " sec" << std::endl; + return 0; +} diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h index 90b128c5eee6..11b7de3ae1b0 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h @@ -68,7 +68,7 @@ class Side_of_oriented_sphere_3 double stz = sz - tz; double st2 = CGAL_NTS square(stx) + CGAL_NTS square(sty) + CGAL_NTS square(stz); - +#if 1 // Compute the semi-static bound. double maxx = CGAL::abs(ptx); double maxy = CGAL::abs(pty); @@ -104,7 +104,11 @@ class Side_of_oriented_sphere_3 if (maxz < artz) maxz = artz; if (maxz < astz) maxz = astz; #endif - +#else +double maxx = (std::max)({CGAL::abs(ptx), CGAL::abs(qtx), CGAL::abs(rtx), CGAL::abs(stx)}); +double maxy = (std::max)({CGAL::abs(pty), CGAL::abs(qty), CGAL::abs(rty), CGAL::abs(sty)}); +double maxz = (std::max)({CGAL::abs(ptz), CGAL::abs(qtz), CGAL::abs(rtz), CGAL::abs(stz)}); +#endif double eps = 1.2466136531027298e-13 * maxx * maxy * maxz; #ifdef CGAL_USE_SSE2_MAX From 7ffb24c5b402b76190f0fb6d518d23a245cdefc5 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 29 Jun 2026 08:48:25 +0100 Subject: [PATCH 2/4] Fix the benchmark progeam so that it does the same number if tests --- .../Filtered_kernel/side_of_oriented_sphere_3.cpp | 12 +++++++----- .../Static_filters/Side_of_oriented_sphere_3.h | 6 ++++++ 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp index 50d83ad74bdc..8dba563251af 100644 --- a/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp +++ b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp @@ -62,14 +62,16 @@ int main(int argc, char* argv[]) std::cout << "We run " << M << " times: All four consecutive points combined with the next " << Q << " points as query" << std::endl; Timer timer; timer.start(); - std::size_t inside = 0; - std::size_t bound = (std::min)(points.size(), std::size_t(Q)); + std::size_t count = 0, inside = 0; + std::size_t bound = (std::min)(points.size()/2, std::size_t(Q)); for(int k = 0; k < M; k++) - for(int i = 0; i < points.size()-6; ++i) - for(int j = i; j < bound ; ++j ) + for(int i = 0; i < points.size()/2; ++i) + for(int j = i+4; j < i+4+bound ; ++j ){ + ++count; if(side_of_oriented_sphere(points[i],points[i+1], points[i+2], points[i+3], points[j]) == CGAL::ON_NEGATIVE_SIDE) ++inside; + } - std ::cout << inside << std::endl; + std ::cout << inside << " inside of " << count << " calls" << std::endl; timer.stop(); std::cout << "Time elapsed: " << timer.time() << " sec" << std::endl; diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h index 11b7de3ae1b0..3be026ee6d4a 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h @@ -108,6 +108,12 @@ class Side_of_oriented_sphere_3 double maxx = (std::max)({CGAL::abs(ptx), CGAL::abs(qtx), CGAL::abs(rtx), CGAL::abs(stx)}); double maxy = (std::max)({CGAL::abs(pty), CGAL::abs(qty), CGAL::abs(rty), CGAL::abs(sty)}); double maxz = (std::max)({CGAL::abs(ptz), CGAL::abs(qtz), CGAL::abs(rtz), CGAL::abs(stz)}); + +/* +double maxx = (std::max)({std::abs(ptx), std::abs(qtx), std::abs(rtx), std::abs(stx)}); +double maxy = (std::max)({std::abs(pty), std::abs(qty), std::abs(rty), std::abs(sty)}); +double maxz = (std::max)({std::abs(ptz), std::abs(qtz), std::abs(rtz), std::abs(stz)}); +*/ #endif double eps = 1.2466136531027298e-13 * maxx * maxy * maxz; From d95499aa1b1b04f133d02e108ff27d2afb47706e Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 29 Jun 2026 12:20:48 +0100 Subject: [PATCH 3/4] Make sure that in the innermost call the 4 points are not always in the same order --- .../benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp index 8dba563251af..1e19c7dcf424 100644 --- a/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp +++ b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp @@ -68,7 +68,7 @@ int main(int argc, char* argv[]) for(int i = 0; i < points.size()/2; ++i) for(int j = i+4; j < i+4+bound ; ++j ){ ++count; - if(side_of_oriented_sphere(points[i],points[i+1], points[i+2], points[i+3], points[j]) == CGAL::ON_NEGATIVE_SIDE) ++inside; + if(side_of_oriented_sphere(points[i+(j%4)],points[i+(1+j%4)], points[i+(2+j%4)], points[i+(3+j%4)], points[j]) == CGAL::ON_NEGATIVE_SIDE) ++inside; } std ::cout << inside << " inside of " << count << " calls" << std::endl; From 05f76c8a916eadcd42ea31c870447014ab9cfa93 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Mon, 29 Jun 2026 21:45:16 +0100 Subject: [PATCH 4/4] Add avx code --- .../benchmark/Filtered_kernel/CMakeLists.txt | 1 + .../side_of_oriented_sphere_3.cpp | 12 ++- .../Side_of_oriented_sphere_3.h | 100 +++++++++++++----- 3 files changed, 84 insertions(+), 29 deletions(-) diff --git a/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt b/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt index f5b0876f25f5..7ae1c7f1549f 100644 --- a/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt +++ b/Filtered_kernel/benchmark/Filtered_kernel/CMakeLists.txt @@ -9,3 +9,4 @@ project(Filtered_kernel_3_Benchmarks) find_package(CGAL REQUIRED) create_single_source_cgal_program("side_of_oriented_sphere_3.cpp") +target_compile_options(side_of_oriented_sphere_3 PRIVATE /arch:AVX2) diff --git a/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp index 1e19c7dcf424..56f672ee18b3 100644 --- a/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp +++ b/Filtered_kernel/benchmark/Filtered_kernel/side_of_oriented_sphere_3.cpp @@ -1,8 +1,12 @@ -// #define CGAL_SMALL_UNORDERED_MAP_STATS // #define CGAL_PROFILE -//#define CGAL_USE_SSE2_FABS -//#define CGAL_USE_SSE2_MAX -//#define CGAL_MSVC_USE_STD_FABS // use this one with precise + +// Nothing defined corresponds to main +// The next two for the three-liner with max(abs*) +// #define CGAL_STD_FABS 1 +// #define CGAL_STD_ABS 1 +// The max(abs*) using AVX +// #define CGAL_VECTORIZE 1 + #define CGAL_NDEBUG 1 #define NDEBUG 1 diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h index 3be026ee6d4a..e2a1dd429ce7 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h @@ -15,9 +15,29 @@ #include #include +#include namespace CGAL { namespace internal { namespace Static_filters_predicates { + + inline double hmax4(__m256d v) +{ + __m128d hi = _mm256_extractf128_pd(v, 1); + __m128d lo = _mm256_castpd256_pd128(v); + __m128d m = _mm_max_pd(lo, hi); + __m128d s = _mm_shuffle_pd(m, m, 1); + m = _mm_max_sd(m, s); + return _mm_cvtsd_f64(m); +} + +inline __m256d abs4(__m256d v) +{ + const __m256d mask = + _mm256_castsi256_pd(_mm256_set1_epi64x(0x7fffffffffffffffLL)); + return _mm256_and_pd(v, mask); +} + + template < typename K_base > class Side_of_oriented_sphere_3 : public K_base::Side_of_oriented_sphere_3 @@ -48,28 +68,68 @@ class Side_of_oriented_sphere_3 { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); +#ifdef CGAL_VECTORIZE + double X[4] = {px-tx, qx - tx, rx - tx, sx - tx}; + double Y[4] = {py-ty, qy - ty, ry - ty, sy - ty}; + double Z[4] = {pz-tz, qz - tz, rz - tz, sz - tz}; + + double pt2 = CGAL_NTS square(X[0]) + CGAL_NTS square(Y[0]) + + CGAL_NTS square(Z[0]); + double qt2 = CGAL_NTS square(X[1]) + CGAL_NTS square(Y[1]) + + CGAL_NTS square(Z[1]); + double rt2 = CGAL_NTS square(X[2]) + CGAL_NTS square(Y[2]) + + CGAL_NTS square(Z[2]); + double st2 = CGAL_NTS square(X[3]) + CGAL_NTS square(Y[3]) + + CGAL_NTS square(Z[3]); + +#else double ptx = px - tx; double pty = py - ty; double ptz = pz - tz; - double pt2 = CGAL_NTS square(ptx) + CGAL_NTS square(pty) - + CGAL_NTS square(ptz); + double qtx = qx - tx; double qty = qy - ty; double qtz = qz - tz; - double qt2 = CGAL_NTS square(qtx) + CGAL_NTS square(qty) - + CGAL_NTS square(qtz); + double rtx = rx - tx; double rty = ry - ty; double rtz = rz - tz; - double rt2 = CGAL_NTS square(rtx) + CGAL_NTS square(rty) - + CGAL_NTS square(rtz); + double stx = sx - tx; double sty = sy - ty; double stz = sz - tz; + + + double pt2 = CGAL_NTS square(ptx) + CGAL_NTS square(pty) + + CGAL_NTS square(ptz); + + double qt2 = CGAL_NTS square(qtx) + CGAL_NTS square(qty) + + CGAL_NTS square(qtz); + double rt2 = CGAL_NTS square(rtx) + CGAL_NTS square(rty) + + CGAL_NTS square(rtz); double st2 = CGAL_NTS square(stx) + CGAL_NTS square(sty) + CGAL_NTS square(stz); -#if 1 +#endif // Compute the semi-static bound. + +#if CGAL_STD_ABS +double maxx = (std::max)({std::abs(ptx), std::abs(qtx), std::abs(rtx), std::abs(stx)}); +double maxy = (std::max)({std::abs(pty), std::abs(qty), std::abs(rty), std::abs(sty)}); +double maxz = (std::max)({std::abs(ptz), std::abs(qtz), std::abs(rtz), std::abs(stz)}); +#elif CGAL_STD_FABS +double maxx = (std::max)({std::fabs(ptx), std::fabs(qtx), std::fabs(rtx), std::fabs(stx)}); +double maxy = (std::max)({std::fabs(pty), std::fabs(qty), std::fabs(rty), std::fabs(sty)}); +double maxz = (std::max)({std::fabs(ptz), std::fabs(qtz), std::fabs(rtz), std::fabs(stz)}); + +#elif CGAL_VECTORIZE +double maxx = hmax4(abs4(_mm256_loadu_pd(X))); +double maxy = hmax4(abs4(_mm256_loadu_pd(Y))); +double maxz = hmax4(abs4(_mm256_loadu_pd(Z))); +#elif CGAL_CGAL_ABS +double maxx = (std::max)({CGAL::abs(ptx), CGAL::abs(qtx), CGAL::abs(rtx), CGAL::abs(stx)}); +double maxy = (std::max)({CGAL::abs(pty), CGAL::abs(qty), CGAL::abs(rty), CGAL::abs(sty)}); +double maxz = (std::max)({CGAL::abs(ptz), CGAL::abs(qtz), CGAL::abs(rtz), CGAL::abs(stz)}); +#else double maxx = CGAL::abs(ptx); double maxy = CGAL::abs(pty); double maxz = CGAL::abs(ptz); @@ -86,12 +146,7 @@ class Side_of_oriented_sphere_3 double artz = CGAL::abs(rtz); double astz = CGAL::abs(stz); -#ifdef CGAL_USE_SSE2_MAX - CGAL::Max mmax; - maxx = mmax(maxx, aqtx, artx, astx); - maxy = mmax(maxy, aqty, arty, asty); - maxz = mmax(maxz, aqtz, artz, astz); -#else + if (maxx < aqtx) maxx = aqtx; if (maxx < artx) maxx = artx; if (maxx < astx) maxx = astx; @@ -103,17 +158,6 @@ class Side_of_oriented_sphere_3 if (maxz < aqtz) maxz = aqtz; if (maxz < artz) maxz = artz; if (maxz < astz) maxz = astz; -#endif -#else -double maxx = (std::max)({CGAL::abs(ptx), CGAL::abs(qtx), CGAL::abs(rtx), CGAL::abs(stx)}); -double maxy = (std::max)({CGAL::abs(pty), CGAL::abs(qty), CGAL::abs(rty), CGAL::abs(sty)}); -double maxz = (std::max)({CGAL::abs(ptz), CGAL::abs(qtz), CGAL::abs(rtz), CGAL::abs(stz)}); - -/* -double maxx = (std::max)({std::abs(ptx), std::abs(qtx), std::abs(rtx), std::abs(stx)}); -double maxy = (std::max)({std::abs(pty), std::abs(qty), std::abs(rty), std::abs(sty)}); -double maxz = (std::max)({std::abs(ptz), std::abs(qtz), std::abs(rtz), std::abs(stz)}); -*/ #endif double eps = 1.2466136531027298e-13 * maxx * maxy * maxz; @@ -136,11 +180,17 @@ double maxz = (std::max)({std::abs(ptz), std::abs(qtz), std::abs(rtz), std::abs( else if (maxy < maxx) std::swap(maxx, maxy); #endif +#ifdef CGAL_VECTORIZE + double det = CGAL::determinant(X[0],Y[0],Z[0],pt2, + X[2],Y[2],Z[2],rt2, + X[1],Y[1],Z[1],qt2, + X[3],Y[3],Z[3],st2); +#else double det = CGAL::determinant(ptx,pty,ptz,pt2, rtx,rty,rtz,rt2, qtx,qty,qtz,qt2, stx,sty,stz,st2); - +#endif // Protect against underflow in the computation of eps. if (maxx < 1e-58) /* sqrt^5(min_double/eps) */ { if (maxx == 0)