1/*
2 Copyright (c) 2005-2019 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 <cstdio>
18#include <vector>
19#include <math.h>
20
21#include "tbb/atomic.h"
22#include "tbb/tick_count.h"
23#include "tbb/task_scheduler_init.h"
24#include "tbb/task_group.h"
25#include "tbb/concurrent_priority_queue.h"
26#include "tbb/spin_mutex.h"
27#include "tbb/parallel_for.h"
28#include "tbb/blocked_range.h"
29#include "../../common/utility/utility.h"
30#include "../../common/utility/fast_random.h"
31
32#if defined(_MSC_VER) && defined(_Wp64)
33 // Workaround for overzealous compiler warnings in /Wp64 mode
34 #pragma warning (disable: 4267)
35#endif /* _MSC_VER && _Wp64 */
36
37using namespace std;
38using namespace tbb;
39
40struct point {
41 double x, y;
42 point() {}
43 point(double _x, double _y) : x(_x), y(_y) {}
44 point(const point& p) : x(p.x), y(p.y) {}
45};
46
47double 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
53point generate_random_point(utility::FastRandom& mr) {
54 const 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
61bool die_toss(size_t a, size_t b, utility::FastRandom& mr) {
62 int node_diff = std::abs((int)(a-b));
63 // near nodes
64 if (node_diff < 16) return true;
65 // mid nodes
66 if (node_diff < 64) return ((int)mr.get() % 8 == 0);
67 // far nodes
68 if (node_diff < 512) return ((int)mr.get() % 16 == 0);
69 return false;
70}
71
72typedef vector<point> point_set;
73typedef size_t vertex_id;
74typedef std::pair<vertex_id,double> vertex_rec;
75typedef vector<vector<vertex_id> > edge_set;
76
77bool verbose = false; // prints bin details and other diagnostics to screen
78bool silent = false; // suppress all output except for time
79size_t N = 1000; // number of vertices
80size_t src = 0; // start of path
81size_t dst = N-1; // end of path
82double INF=100000.0; // infinity
83size_t grainsize = 16; // number of vertices per task on average
84size_t max_spawn; // max tasks to spawn
85tbb::atomic<size_t> num_spawn; // number of active tasks
86
87point_set vertices; // vertices
88edge_set edges; // edges
89vector<vertex_id> predecessor; // for recreating path from src to dst
90
91vector<double> f_distance; // estimated distances at particular vertex
92vector<double> g_distance; // current shortest distances from src vertex
93spin_mutex *locks; // a lock for each vertex
94task_group *sp_group; // task group for tasks executing sub-problems
95
96class compare_f {
97public:
98 bool operator()(const vertex_rec& u, const vertex_rec& v) const {
99 return u.second>v.second;
100 }
101};
102
103concurrent_priority_queue<vertex_rec, compare_f> open_set; // tentative vertices
104
105void shortpath_helper();
106
107#if !__TBB_CPP11_LAMBDAS_PRESENT
108class shortpath_helper_functor {
109public:
110 shortpath_helper_functor() {};
111 void operator() () const { shortpath_helper(); }
112};
113#endif
114
115void shortpath() {
116 sp_group = new task_group;
117 g_distance[src] = 0.0; // src's distance from src is zero
118 f_distance[src] = get_distance(vertices[src], vertices[dst]); // estimate distance from src to dst
119 open_set.push(make_pair(src,f_distance[src])); // push src into open_set
120#if __TBB_CPP11_LAMBDAS_PRESENT
121 sp_group->run([](){ shortpath_helper(); });
122#else
123 sp_group->run( shortpath_helper_functor() );
124#endif
125 sp_group->wait();
126 delete sp_group;
127}
128
129void shortpath_helper() {
130 vertex_rec u_rec;
131 while (open_set.try_pop(u_rec)) {
132 vertex_id u = u_rec.first;
133 if (u==dst) continue;
134 double f = u_rec.second;
135 double old_g_u = 0.0;
136 {
137 spin_mutex::scoped_lock l(locks[u]);
138 if (f > f_distance[u]) continue; // prune search space
139 old_g_u = g_distance[u];
140 }
141 for (size_t i=0; i<edges[u].size(); ++i) {
142 vertex_id v = edges[u][i];
143 double new_g_v = old_g_u + get_distance(vertices[u], vertices[v]);
144 double new_f_v = 0.0;
145 // the push flag lets us move some work out of the critical section below
146 bool push = false;
147 {
148 spin_mutex::scoped_lock l(locks[v]);
149 if (new_g_v < g_distance[v]) {
150 predecessor[v] = u;
151 g_distance[v] = new_g_v;
152 new_f_v = f_distance[v] = g_distance[v] + get_distance(vertices[v], vertices[dst]);
153 push = true;
154 }
155 }
156 if (push) {
157 open_set.push(make_pair(v,new_f_v));
158 size_t n_spawn = ++num_spawn;
159 if (n_spawn < max_spawn) {
160#if __TBB_CPP11_LAMBDAS_PRESENT
161 sp_group->run([]{ shortpath_helper(); });
162#else
163 sp_group->run( shortpath_helper_functor() );
164#endif
165 }
166 else --num_spawn;
167 }
168 }
169 }
170 --num_spawn;
171}
172
173void make_path(vertex_id src, vertex_id dst, vector<vertex_id>& path) {
174 vertex_id at = predecessor[dst];
175 if (at == N) path.push_back(src);
176 else if (at == src) { path.push_back(src); path.push_back(dst); }
177 else { make_path(src, at, path); path.push_back(dst); }
178}
179
180void print_path() {
181 vector<vertex_id> path;
182 double path_length=0.0;
183 make_path(src, dst, path);
184 if (verbose) printf("\n ");
185 for (size_t i=0; i<path.size(); ++i) {
186 if (path[i] != dst) {
187 double seg_length = get_distance(vertices[path[i]], vertices[path[i+1]]);
188 if (verbose) printf("%6.1f ", seg_length);
189 path_length += seg_length;
190 }
191 else if (verbose) printf("\n");
192 }
193 if (verbose) {
194 for (size_t i=0; i<path.size(); ++i) {
195 if (path[i] != dst) printf("(%4d)------>", (int)path[i]);
196 else printf("(%4d)\n", (int)path[i]);
197 }
198 }
199 if (verbose) printf("Total distance = %5.1f\n", path_length);
200 else if (!silent) printf(" %5.1f\n", path_length);
201}
202
203int get_default_num_threads() {
204 static int threads = 0;
205 if (threads == 0)
206 threads = tbb::task_scheduler_init::default_num_threads();
207 return threads;
208}
209
210#if !__TBB_CPP11_LAMBDAS_PRESENT
211class gen_vertices {
212public:
213 gen_vertices() {}
214 void operator() (blocked_range<size_t>& r) const {
215 utility::FastRandom my_random((unsigned int)r.begin());
216 for (size_t i=r.begin(); i!=r.end(); ++i) {
217 vertices[i] = generate_random_point(my_random);
218 }
219 }
220};
221
222class gen_edges {
223public:
224 gen_edges() {}
225 void operator() (blocked_range<size_t>& r) const {
226 utility::FastRandom my_random((unsigned int)r.begin());
227 for (size_t i=r.begin(); i!=r.end(); ++i) {
228 for (size_t j=0; j<i; ++j) {
229 if (die_toss(i, j, my_random))
230 edges[i].push_back(j);
231 }
232 }
233 }
234};
235
236class reset_vertices {
237public:
238 reset_vertices() {}
239 void operator() (blocked_range<size_t>& r) const {
240 for (size_t i=r.begin(); i!=r.end(); ++i) {
241 f_distance[i] = g_distance[i] = INF;
242 predecessor[i] = N;
243 }
244 }
245};
246#endif
247
248void InitializeGraph() {
249 task_scheduler_init init(get_default_num_threads());
250 vertices.resize(N);
251 edges.resize(N);
252 predecessor.resize(N);
253 g_distance.resize(N);
254 f_distance.resize(N);
255 locks = new spin_mutex[N];
256 if (verbose) printf("Generating vertices...\n");
257#if __TBB_CPP11_LAMBDAS_PRESENT
258 parallel_for(blocked_range<size_t>(0,N,64),
259 [&](blocked_range<size_t>& r) {
260 utility::FastRandom my_random(r.begin());
261 for (size_t i=r.begin(); i!=r.end(); ++i) {
262 vertices[i] = generate_random_point(my_random);
263 }
264 }, simple_partitioner());
265#else
266 parallel_for(blocked_range<size_t>(0,N,64), gen_vertices(), simple_partitioner());
267#endif
268 if (verbose) printf("Generating edges...\n");
269#if __TBB_CPP11_LAMBDAS_PRESENT
270 parallel_for(blocked_range<size_t>(0,N,64),
271 [&](blocked_range<size_t>& r) {
272 utility::FastRandom my_random(r.begin());
273 for (size_t i=r.begin(); i!=r.end(); ++i) {
274 for (size_t j=0; j<i; ++j) {
275 if (die_toss(i, j, my_random))
276 edges[i].push_back(j);
277 }
278 }
279 }, simple_partitioner());
280#else
281 parallel_for(blocked_range<size_t>(0,N,64), gen_edges(), simple_partitioner());
282#endif
283 for (size_t i=0; i<N; ++i) {
284 for (size_t j=0; j<edges[i].size(); ++j) {
285 vertex_id k = edges[i][j];
286 edges[k].push_back(i);
287 }
288 }
289 if (verbose) printf("Done.\n");
290}
291
292void ReleaseGraph() {
293 delete []locks;
294}
295
296void ResetGraph() {
297 task_scheduler_init init(get_default_num_threads());
298#if __TBB_CPP11_LAMBDAS_PRESENT
299 parallel_for(blocked_range<size_t>(0,N),
300 [&](blocked_range<size_t>& r) {
301 for (size_t i=r.begin(); i!=r.end(); ++i) {
302 f_distance[i] = g_distance[i] = INF;
303 predecessor[i] = N;
304 }
305 });
306#else
307 parallel_for(blocked_range<size_t>(0,N), reset_vertices());
308#endif
309}
310
311int main(int argc, char *argv[]) {
312 try {
313 utility::thread_number_range threads(get_default_num_threads);
314 utility::parse_cli_arguments(argc, argv,
315 utility::cli_argument_pack()
316 //"-h" option for displaying help is present implicitly
317 .positional_arg(threads,"#threads",utility::thread_number_range_desc)
318 .arg(verbose,"verbose"," print diagnostic output to screen")
319 .arg(silent,"silent"," limits output to timing info; overrides verbose")
320 .arg(N,"N"," number of vertices")
321 .arg(src,"start"," start of path")
322 .arg(dst,"end"," end of path")
323 );
324 if (silent) verbose = false; // make silent override verbose
325 else
326 printf("shortpath will run with %d vertices to find shortest path between vertices"
327 " %d and %d using %d:%d threads.\n",
328 (int)N, (int)src, (int)dst, (int)threads.first, (int)threads.last);
329
330 if (dst >= N) {
331 if (verbose)
332 printf("end value %d is invalid for %d vertices; correcting to %d\n", (int)dst, (int)N, (int)N-1);
333 dst = N-1;
334 }
335
336 num_spawn = 0;
337 max_spawn = N/grainsize;
338 tick_count t0, t1;
339 InitializeGraph();
340 for (int n_thr=threads.first; n_thr<=threads.last; n_thr=threads.step(n_thr)) {
341 ResetGraph();
342 task_scheduler_init init(n_thr);
343 t0 = tick_count::now();
344 shortpath();
345 t1 = tick_count::now();
346 if (!silent) {
347 if (predecessor[dst] != N) {
348 printf("%d threads: [%6.6f] The shortest path from vertex %d to vertex %d is:",
349 (int)n_thr, (t1-t0).seconds(), (int)src, (int)dst);
350 print_path();
351 }
352 else {
353 printf("%d threads: [%6.6f] There is no path from vertex %d to vertex %d\n",
354 (int)n_thr, (t1-t0).seconds(), (int)src, (int)dst);
355 }
356 } else
357 utility::report_elapsed_time((t1-t0).seconds());
358 }
359 ReleaseGraph();
360 return 0;
361 } catch(std::exception& e) {
362 cerr<<"error occurred. error text is :\"" <<e.what()<<"\"\n";
363 return 1;
364 }
365}
366