Read-safe multithreading with Lazy_exact_nt #2685
Description
I can't seem to find any discussion in the documentation the thread safety (or lack thereof) when using CGAL:: Lazy_exact_nt, e.g., via CGAL::Epeck. I'm worried that this type uses non-atomic reference counting implying that even reading a variable could result in a race condition and subsequent memory violation.
In my specific case, I have a giant list of points that should be reduced to constrained Delaunay triangulations (in my simplified code below, I put every 4 points into a CDT). Conceptually, my problem is embarassingly parallel with respect to the "points". However, my points are constructed via arithmetic involving common points, so I understand that the lazy evaluation makes this a non-trivial dependency graph (simulated below by a simple addition of the first point to all others). Still, my problem is completely parallel with respect to the outputs: no CDT depends on any other one.
Here's a minimal-ish reproducible example:
#define SERIAL
#include <vector>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Vector_2.h>
#include <CGAL/Point_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Constrained_triangulation_plus_2.h>
#include <cmath>
#include <cassert>
#include <thread>
#include <vector>
#include <algorithm>
template<
typename Index,
typename FunctionType>
inline bool parallel_for(
const Index loop_size,
const FunctionType & func)
{
assert(loop_size>=0);
if(loop_size==0) return false;
// Estimate number of threads in the pool
// http://ideone.com/Z7zldb
const static size_t sthc = std::thread::hardware_concurrency();
const size_t nthreads = (sthc==0?8:sthc);
// Size of a slice for the range functions
Index slice =
std::max(
(Index)std::round((loop_size+1)/static_cast<double>(nthreads)),(Index)1);
// [Helper] Inner loop
const auto & range = [&func](const Index k1, const Index k2, const size_t t)
{
for(Index k = k1; k < k2; k++) func(k);
};
// Create pool and launch jobs
std::vector<std::thread> pool;
pool.reserve(nthreads);
// Inner range extents
Index i1 = 0;
Index i2 = std::min(0 + slice, loop_size);
{
size_t t = 0;
for (; t+1 < nthreads && i1 < loop_size; ++t)
{
pool.emplace_back(range, i1, i2, t);
i1 = i2;
i2 = std::min(i2 + slice, loop_size);
}
if (i1 < loop_size)
{
pool.emplace_back(range, i1, loop_size, t);
}
}
// Wait for jobs to finish
for (std::thread &t : pool) if (t.joinable()) t.join();
return true;
}
int main(int argc, char *argv[])
{
typedef CGAL::Epeck Kernel;
typedef CGAL::Triangulation_vertex_base_2<Kernel> TVB_2;
typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
typedef CGAL::Exact_intersections_tag Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> CDT_2;
typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
std::vector< CGAL::Point_2<CGAL::Epeck> > points(100000);
points[0] = CGAL::Point_2<CGAL::Epeck>(0,1);
for(int i = 1;i<points.size();i++)
{
points[i] = points[(i-1)%5]+(CGAL::Vector_2<CGAL::Epeck>(i%7,i%11))/CGAL::Epeck::FT(3);
}
std::vector< CDT_plus_2 > cdts(points.size()/4);
#ifdef SERIAL
for(int i = 1;i<cdts.size();i++)
#else
parallel_for(cdts.size(),[&](int i)
#endif
{
CDT_plus_2 cdt;
cdt.insert_constraint(points[4*i+0],points[4*i+1]);
cdt.insert_constraint(points[4*i+2],points[4*i+3]);
}
#ifndef SERIAL
);
#endif
}Use #define SERIAL to control whether the code executes in serial (should exit cleanly) or in parallel (usually causes EXC_BAD_ACCESS fault).
FWIW @qnzhou expressed a thought that perhaps this reference counting is easy to fix using atomic integers
Environment
- Operating system (Windows/Mac/Linux, 32/64 bits): Mac
- Compiler: Apple LLVM version 9.0.0 (clang-900.0.39.2)
- Release or debug mode: either
- Specific flags used (if any):
- CGAL version: 1040911000
- Boost version: 106400
- Other libraries versions if used (Eigen, TBB, etc.): (std::thread-based parallel_for borrowed from libigl)