Skip to content

Commit 534daa9

Browse files
zferguszfergus
committed
Replace scalar enumerable_thread_specific with parallel_reduce
1 parent 86f1db3 commit 534daa9

File tree

4 files changed

+48
-57
lines changed

4 files changed

+48
-57
lines changed

‎src/ipc/collisions/collisions.cpp

+10-13
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include <ipc/utils/local_to_global.hpp>
1010

1111
#include <tbb/parallel_for.h>
12+
#include <tbb/parallel_reduce.h>
1213
#include <tbb/parallel_sort.h>
1314
#include <tbb/blocked_range.h>
1415
#include <tbb/enumerable_thread_specific.h>
@@ -328,25 +329,21 @@ double Collisions::compute_minimum_distance(
328329
const Eigen::MatrixXi& edges = mesh.edges();
329330
const Eigen::MatrixXi& faces = mesh.faces();
330331

331-
tbb::enumerable_thread_specific<double> storage(
332-
std::numeric_limits<double>::infinity());
333-
334-
tbb::parallel_for(
335-
tbb::blocked_range<size_t>(0, size()),
336-
[&](tbb::blocked_range<size_t> r) {
337-
double& local_min_dist = storage.local();
338-
332+
return tbb::parallel_reduce(
333+
/*range=*/tbb::blocked_range<size_t>(0, size()),
334+
/*initial_min_distance=*/std::numeric_limits<double>::infinity(),
335+
[&](tbb::blocked_range<size_t> r, double partial_min_dist) -> double {
339336
for (size_t i = r.begin(); i < r.end(); i++) {
340337
const double dist = (*this)[i].compute_distance(
341338
(*this)[i].dof(vertices, edges, faces));
342339

343-
if (dist < local_min_dist) {
344-
local_min_dist = dist;
340+
if (dist < partial_min_dist) {
341+
partial_min_dist = dist;
345342
}
346343
}
347-
});
348-
349-
return storage.combine([](double a, double b) { return std::min(a, b); });
344+
return partial_min_dist;
345+
},
346+
/*combine=*/[](double a, double b) { return std::min(a, b); });
350347
}
351348

352349
// ============================================================================

‎src/ipc/implicits/plane.cpp

+7-11
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,8 @@
33
#include <ipc/distance/point_plane.hpp>
44
#include <ipc/ccd/point_static_plane.hpp>
55

6-
#include <tbb/parallel_for.h>
6+
#include <tbb/parallel_reduce.h>
77
#include <tbb/blocked_range.h>
8-
#include <tbb/enumerable_thread_specific.h>
98

109
namespace ipc {
1110

@@ -98,13 +97,10 @@ double compute_point_plane_collision_free_stepsize(
9897
assert(plane_normals.rows() == n_planes);
9998
assert(points_t0.rows() == points_t1.rows());
10099

101-
tbb::enumerable_thread_specific<double> storage(1);
102-
103-
tbb::parallel_for(
100+
const double earliest_toi = tbb::parallel_reduce(
104101
tbb::blocked_range<size_t>(0, points_t0.rows()),
105-
[&](tbb::blocked_range<size_t> r) {
106-
double& earliest_toi = storage.local();
107-
102+
/*inital_step_size=*/1.0,
103+
[&](tbb::blocked_range<size_t> r, double earliest_toi) {
108104
for (size_t vi = r.begin(); vi < r.end(); vi++) {
109105
for (size_t pi = 0; pi < n_planes; pi++) {
110106
if (!can_collide(vi, pi)) {
@@ -126,10 +122,10 @@ double compute_point_plane_collision_free_stepsize(
126122
}
127123
}
128124
}
129-
});
125+
return earliest_toi;
126+
},
127+
[&](double a, double b) { return std::min(a, b); });
130128

131-
const double earliest_toi =
132-
storage.combine([](double a, double b) { return std::min(a, b); });
133129
assert(earliest_toi >= 0 && earliest_toi <= 1.0);
134130
return earliest_toi;
135131
}

‎src/ipc/potentials/friction_potential.cpp

+2-3
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ Eigen::VectorXd FrictionPotential::force(
2828
const Eigen::MatrixXi& edges = mesh.edges();
2929
const Eigen::MatrixXi& faces = mesh.faces();
3030

31-
tbb::enumerable_thread_specific<Eigen::VectorXd> storage(
31+
tbb::combinable<Eigen::VectorXd> storage(
3232
Eigen::VectorXd::Zero(velocities.size()));
3333

3434
tbb::parallel_for(
@@ -52,8 +52,7 @@ Eigen::VectorXd FrictionPotential::force(
5252
}
5353
});
5454

55-
return storage.combine([](const Eigen::VectorXd& a,
56-
const Eigen::VectorXd& b) { return a + b; });
55+
return storage.combine(std::plus<const Eigen::VectorXd&>());
5756
}
5857

5958
Eigen::SparseMatrix<double> FrictionPotential::force_jacobian(

‎src/ipc/potentials/potential.tpp

+29-30
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,11 @@
44

55
#include <ipc/utils/local_to_global.hpp>
66

7-
#include <tbb/parallel_for.h>
87
#include <tbb/blocked_range.h>
8+
#include <tbb/combinable.h>
99
#include <tbb/enumerable_thread_specific.h>
10+
#include <tbb/parallel_for.h>
11+
#include <tbb/parallel_reduce.h>
1012

1113
namespace ipc {
1214

@@ -18,25 +20,18 @@ double Potential<TCollisions>::operator()(
1820
{
1921
assert(X.rows() == mesh.num_vertices());
2022

21-
if (collisions.empty()) {
22-
return 0;
23-
}
24-
25-
tbb::enumerable_thread_specific<double> storage(0);
26-
27-
tbb::parallel_for(
28-
tbb::blocked_range<size_t>(size_t(0), collisions.size()),
29-
[&](const tbb::blocked_range<size_t>& r) {
30-
auto& local_potential = storage.local();
23+
return tbb::parallel_reduce(
24+
tbb::blocked_range<size_t>(size_t(0), collisions.size()), 0.0,
25+
[&](const tbb::blocked_range<size_t>& r, double partial_sum) {
3126
for (size_t i = r.begin(); i < r.end(); i++) {
3227
// Quadrature weight is premultiplied by local potential
33-
local_potential += (*this)(
28+
partial_sum += (*this)(
3429
collisions[i],
3530
collisions[i].dof(X, mesh.edges(), mesh.faces()));
3631
}
37-
});
38-
39-
return storage.combine([](double a, double b) { return a + b; });
32+
return partial_sum;
33+
},
34+
std::plus<double>());
4035
}
4136

4237
template <class TCollisions>
@@ -53,14 +48,11 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
5348

5449
const int dim = X.cols();
5550

56-
tbb::enumerable_thread_specific<Eigen::VectorXd> storage(
57-
Eigen::VectorXd::Zero(X.size()));
51+
tbb::combinable<Eigen::VectorXd> grad(Eigen::VectorXd::Zero(X.size()));
5852

5953
tbb::parallel_for(
6054
tbb::blocked_range<size_t>(size_t(0), collisions.size()),
6155
[&](const tbb::blocked_range<size_t>& r) {
62-
auto& global_grad = storage.local();
63-
6456
for (size_t i = r.begin(); i < r.end(); i++) {
6557
const TCollision& collision = collisions[i];
6658

@@ -71,12 +63,11 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
7163
collision.vertex_ids(mesh.edges(), mesh.faces());
7264

7365
local_gradient_to_global_gradient(
74-
local_grad, vids, dim, global_grad);
66+
local_grad, vids, dim, grad.local());
7567
}
7668
});
7769

78-
return storage.combine([](const Eigen::VectorXd& a,
79-
const Eigen::VectorXd& b) { return a + b; });
70+
return grad.combine(std::plus<const Eigen::VectorXd&>());
8071
}
8172

8273
template <class TCollisions>
@@ -121,14 +112,22 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
121112
}
122113
});
123114

124-
Eigen::SparseMatrix<double> hess(ndof, ndof);
125-
for (const auto& local_hess_triplets : storage) {
126-
Eigen::SparseMatrix<double> local_hess(ndof, ndof);
127-
local_hess.setFromTriplets(
128-
local_hess_triplets.begin(), local_hess_triplets.end());
129-
hess += local_hess;
130-
}
131-
return hess;
115+
// Combine the local hessians
116+
tbb::combinable<Eigen::SparseMatrix<double>> hess(
117+
Eigen::SparseMatrix<double>(ndof, ndof));
118+
119+
tbb::parallel_for(
120+
tbb::blocked_range<decltype(storage)::iterator>(
121+
storage.begin(), storage.end()),
122+
[&](const tbb::blocked_range<decltype(storage)::iterator>& r) {
123+
for (auto it = r.begin(); it != r.end(); ++it) {
124+
Eigen::SparseMatrix<double> local_hess(ndof, ndof);
125+
local_hess.setFromTriplets(it->begin(), it->end());
126+
hess.local() += local_hess;
127+
}
128+
});
129+
130+
return hess.combine(std::plus<const Eigen::SparseMatrix<double>&>());
132131
}
133132

134133
} // namespace ipc

0 commit comments

Comments
 (0)