1 /*
2 Copyright (c) 2005-2021 Intel Corporation
3
4 Licensed under the Apache License, Version 2.0 (the "License");
5 you may not use this file except in compliance with the License.
6 You may obtain a copy of the License at
7
8 http://www.apache.org/licenses/LICENSE-2.0
9
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
15 */
16
17 #include <cmath>
18 #include <cstdio>
19
20 #include <vector>
21 #include <atomic>
22
23 #include "oneapi/tbb/tick_count.h"
24 #include "oneapi/tbb/task_group.h"
25 #include "oneapi/tbb/concurrent_priority_queue.h"
26 #include "oneapi/tbb/spin_mutex.h"
27 #include "oneapi/tbb/parallel_for.h"
28 #include "oneapi/tbb/blocked_range.h"
29 #include "oneapi/tbb/global_control.h"
30
31 #include "common/utility/utility.hpp"
32 #include "common/utility/fast_random.hpp"
33 #include "common/utility/get_default_num_threads.hpp"
34
35 #if defined(_MSC_VER) && defined(_Wp64)
36 // Workaround for overzealous compiler warnings in /Wp64 mode
37 #pragma warning(disable : 4267)
38 #endif /* _MSC_VER && _Wp64 */
39
40 struct point {
41 double x, y;
pointpoint42 point() {}
pointpoint43 point(double _x, double _y) : x(_x), y(_y) {}
pointpoint44 point(const point& p) : x(p.x), y(p.y) {}
45 };
46
get_distance(const point & p1,const point & p2)47 double get_distance(const point& p1, const point& p2) {
48 double xdiff = p1.x - p2.x, ydiff = p1.y - p2.y;
49 return sqrt(xdiff * xdiff + ydiff * ydiff);
50 }
51
52 // generates random points on 2D plane within a box of maxsize width & height
generate_random_point(utility::FastRandom & mr)53 point generate_random_point(utility::FastRandom& mr) {
54 const std::size_t maxsize = 500;
55 double x = (double)(mr.get() % maxsize);
56 double y = (double)(mr.get() % maxsize);
57 return point(x, y);
58 }
59
60 // weighted toss makes closer nodes (in the point vector) heavily connected
die_toss(std::size_t a,std::size_t b,utility::FastRandom & mr)61 bool die_toss(std::size_t a, std::size_t b, utility::FastRandom& mr) {
62 int node_diff = std::abs(int(a - b));
63 // near nodes
64 if (node_diff < 16)
65 return true;
66 // mid nodes
67 if (node_diff < 64)
68 return ((int)mr.get() % 8 == 0);
69 // far nodes
70 if (node_diff < 512)
71 return ((int)mr.get() % 16 == 0);
72 return false;
73 }
74
75 typedef std::vector<point> point_set;
76 typedef std::size_t vertex_id;
77 typedef std::pair<vertex_id, double> vertex_rec;
78 typedef std::vector<std::vector<vertex_id>> edge_set;
79
80 bool verbose = false; // prints bin details and other diagnostics to screen
81 bool silent = false; // suppress all output except for time
82 std::size_t N = 1000; // number of vertices
83 std::size_t src = 0; // start of path
84 std::size_t dst = N - 1; // end of path
85 double INF = 100000.0; // infinity
86 std::size_t grainsize = 16; // number of vertices per task on average
87 std::size_t max_spawn; // max tasks to spawn
88 std::atomic<std::size_t> num_spawn; // number of active tasks
89
90 point_set vertices; // vertices
91 edge_set edges; // edges
92 std::vector<vertex_id> predecessor; // for recreating path from src to dst
93
94 std::vector<double> f_distance; // estimated distances at particular vertex
95 std::vector<double> g_distance; // current shortest distances from src vertex
96 oneapi::tbb::spin_mutex* locks; // a lock for each vertex
97 oneapi::tbb::task_group* sp_group; // task group for tasks executing sub-problems
98
99 struct compare_f {
operator ()compare_f100 bool operator()(const vertex_rec& u, const vertex_rec& v) const {
101 return u.second > v.second;
102 }
103 };
104
105 oneapi::tbb::concurrent_priority_queue<vertex_rec, compare_f> open_set; // tentative vertices
106
107 void shortpath_helper();
108
shortpath()109 void shortpath() {
110 sp_group = new oneapi::tbb::task_group;
111 g_distance[src] = 0.0; // src's distance from src is zero
112 f_distance[src] =
113 get_distance(vertices[src], vertices[dst]); // estimate distance from src to dst
114 open_set.emplace(src, f_distance[src]); // emplace src into open_set
115 sp_group->run([]() {
116 shortpath_helper();
117 });
118 sp_group->wait();
119 delete sp_group;
120 }
121
shortpath_helper()122 void shortpath_helper() {
123 vertex_rec u_rec;
124 while (open_set.try_pop(u_rec)) {
125 vertex_id u = u_rec.first;
126 if (u == dst)
127 continue;
128 double f = u_rec.second;
129 double old_g_u = 0.0;
130 {
131 oneapi::tbb::spin_mutex::scoped_lock l(locks[u]);
132 if (f > f_distance[u])
133 continue; // prune search space
134 old_g_u = g_distance[u];
135 }
136 for (std::size_t i = 0; i < edges[u].size(); ++i) {
137 vertex_id v = edges[u][i];
138 double new_g_v = old_g_u + get_distance(vertices[u], vertices[v]);
139 double new_f_v = 0.0;
140 // the push flag lets us move some work out of the critical section below
141 bool push = false;
142 {
143 oneapi::tbb::spin_mutex::scoped_lock l(locks[v]);
144 if (new_g_v < g_distance[v]) {
145 predecessor[v] = u;
146 g_distance[v] = new_g_v;
147 new_f_v = f_distance[v] =
148 g_distance[v] + get_distance(vertices[v], vertices[dst]);
149 push = true;
150 }
151 }
152 if (push) {
153 open_set.push(std::make_pair(v, new_f_v));
154 std::size_t n_spawn = ++num_spawn;
155 if (n_spawn < max_spawn) {
156 sp_group->run([] {
157 shortpath_helper();
158 });
159 }
160 else
161 --num_spawn;
162 }
163 }
164 }
165 --num_spawn;
166 }
167
make_path(vertex_id src,vertex_id dst,std::vector<vertex_id> & path)168 void make_path(vertex_id src, vertex_id dst, std::vector<vertex_id>& path) {
169 vertex_id at = predecessor[dst];
170 if (at == N)
171 path.push_back(src);
172 else if (at == src) {
173 path.push_back(src);
174 path.push_back(dst);
175 }
176 else {
177 make_path(src, at, path);
178 path.push_back(dst);
179 }
180 }
181
print_path()182 void print_path() {
183 std::vector<vertex_id> path;
184 double path_length = 0.0;
185 make_path(src, dst, path);
186 if (verbose)
187 printf("\n ");
188 for (std::size_t i = 0; i < path.size(); ++i) {
189 if (path[i] != dst) {
190 double seg_length = get_distance(vertices[path[i]], vertices[path[i + 1]]);
191 if (verbose)
192 printf("%6.1f ", seg_length);
193 path_length += seg_length;
194 }
195 else if (verbose)
196 printf("\n");
197 }
198 if (verbose) {
199 for (std::size_t i = 0; i < path.size(); ++i) {
200 if (path[i] != dst)
201 printf("(%4d)------>", (int)path[i]);
202 else
203 printf("(%4d)\n", (int)path[i]);
204 }
205 }
206 if (verbose)
207 printf("Total distance = %5.1f\n", path_length);
208 else if (!silent)
209 printf(" %5.1f\n", path_length);
210 }
211
InitializeGraph()212 void InitializeGraph() {
213 oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism,
214 utility::get_default_num_threads());
215 vertices.resize(N);
216 edges.resize(N);
217 predecessor.resize(N);
218 g_distance.resize(N);
219 f_distance.resize(N);
220 locks = new oneapi::tbb::spin_mutex[N];
221
222 if (verbose)
223 printf("Generating vertices...\n");
224 oneapi::tbb::parallel_for(
225 oneapi::tbb::blocked_range<std::size_t>(0, N, 64),
226 [&](oneapi::tbb::blocked_range<std::size_t>& r) {
227 utility::FastRandom my_random(r.begin());
228 for (std::size_t i = r.begin(); i != r.end(); ++i) {
229 vertices[i] = generate_random_point(my_random);
230 }
231 },
232 oneapi::tbb::simple_partitioner());
233
234 if (verbose)
235 printf("Generating edges...\n");
236 oneapi::tbb::parallel_for(
237 oneapi::tbb::blocked_range<std::size_t>(0, N, 64),
238 [&](oneapi::tbb::blocked_range<std::size_t>& r) {
239 utility::FastRandom my_random(r.begin());
240 for (std::size_t i = r.begin(); i != r.end(); ++i) {
241 for (std::size_t j = 0; j < i; ++j) {
242 if (die_toss(i, j, my_random))
243 edges[i].push_back(j);
244 }
245 }
246 },
247 oneapi::tbb::simple_partitioner());
248
249 for (std::size_t i = 0; i < N; ++i) {
250 for (std::size_t j = 0; j < edges[i].size(); ++j) {
251 vertex_id k = edges[i][j];
252 edges[k].push_back(i);
253 }
254 }
255 if (verbose)
256 printf("Done.\n");
257 }
258
ReleaseGraph()259 void ReleaseGraph() {
260 delete[] locks;
261 }
262
ResetGraph()263 void ResetGraph() {
264 oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism,
265 utility::get_default_num_threads());
266 oneapi::tbb::parallel_for(oneapi::tbb::blocked_range<std::size_t>(0, N),
267 [&](oneapi::tbb::blocked_range<std::size_t>& r) {
268 for (std::size_t i = r.begin(); i != r.end(); ++i) {
269 f_distance[i] = g_distance[i] = INF;
270 predecessor[i] = N;
271 }
272 });
273 }
274
main(int argc,char * argv[])275 int main(int argc, char* argv[]) {
276 utility::thread_number_range threads(utility::get_default_num_threads);
277 utility::parse_cli_arguments(
278 argc,
279 argv,
280 utility::cli_argument_pack()
281 //"-h" option for displaying help is present implicitly
282 .positional_arg(threads, "#threads", utility::thread_number_range_desc)
283 .arg(verbose, "verbose", " print diagnostic output to screen")
284 .arg(silent, "silent", " limits output to timing info; overrides verbose")
285 .arg(N, "N", " number of vertices")
286 .arg(src, "start", " start of path")
287 .arg(dst, "end", " end of path"));
288 if (silent)
289 verbose = false; // make silent override verbose
290 else
291 printf("shortpath will run with %d vertices to find shortest path between vertices"
292 " %d and %d using %d:%d threads.\n",
293 (int)N,
294 (int)src,
295 (int)dst,
296 (int)threads.first,
297 (int)threads.last);
298
299 if (dst >= N) {
300 if (verbose)
301 printf("end value %d is invalid for %d vertices; correcting to %d\n",
302 (int)dst,
303 (int)N,
304 (int)N - 1);
305 dst = N - 1;
306 }
307
308 num_spawn = 0;
309 max_spawn = N / grainsize;
310 oneapi::tbb::tick_count t0, t1;
311 InitializeGraph();
312 for (int n_thr = threads.first; n_thr <= threads.last; n_thr = threads.step(n_thr)) {
313 ResetGraph();
314 oneapi::tbb::global_control c(oneapi::tbb::global_control::max_allowed_parallelism, n_thr);
315 t0 = oneapi::tbb::tick_count::now();
316 shortpath();
317 t1 = oneapi::tbb::tick_count::now();
318 if (!silent) {
319 if (predecessor[dst] != N) {
320 printf("%d threads: [%6.6f] The shortest path from vertex %d to vertex %d is:",
321 (int)n_thr,
322 (t1 - t0).seconds(),
323 (int)src,
324 (int)dst);
325 print_path();
326 }
327 else {
328 printf("%d threads: [%6.6f] There is no path from vertex %d to vertex %d\n",
329 (int)n_thr,
330 (t1 - t0).seconds(),
331 (int)src,
332 (int)dst);
333 }
334 }
335 else
336 utility::report_elapsed_time((t1 - t0).seconds());
337 }
338 ReleaseGraph();
339 return 0;
340 }
341