examples/concurrent_priority_queue/shortpath/shortpath.cpp

   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
  37 using namespace std;
  38 using namespace tbb;
  39
  40 struct 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
  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
  53 point 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
  61 bool 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
  72 typedef vector<point> point_set;
  73 typedef size_t vertex_id;
  74 typedef std::pair<vertex_id,double> vertex_rec;
  75 typedef vector<vector<vertex_id> > edge_set;
  76
  77 bool verbose = false;          // prints bin details and other diagnostics to screen
  78 bool silent = false;           // suppress all output except for time
  79 size_t N = 1000;               // number of vertices
  80 size_t src = 0;                // start of path
  81 size_t dst = N-1;              // end of path
  82 double INF=100000.0;           // infinity
  83 size_t grainsize = 16;         // number of vertices per task on average
  84 size_t max_spawn;              // max tasks to spawn
  85 tbb::atomic<size_t> num_spawn;      // number of active tasks
  86
  87 point_set vertices;            // vertices
  88 edge_set edges;                // edges
  89 vector<vertex_id> predecessor; // for recreating path from src to dst
  90
  91 vector<double> f_distance;     // estimated distances at particular vertex
  92 vector<double> g_distance;     // current shortest distances from src vertex
  93 spin_mutex    *locks;          // a lock for each vertex
  94 task_group *sp_group;          // task group for tasks executing sub-problems
  95
  96 class compare_f {
  97 public:
  98     bool operator()(const vertex_rec& u, const vertex_rec& v) const {
  99         return u.second>v.second;
 100     }
 101 };
 102
 103 concurrent_priority_queue<vertex_rec, compare_f> open_set; // tentative vertices
 104
 105 void shortpath_helper();
 106
 107 #if !__TBB_CPP11_LAMBDAS_PRESENT
 108 class shortpath_helper_functor {
 109 public:
 110     shortpath_helper_functor() {};
 111     void operator() () const { shortpath_helper(); }
 112 };
 113 #endif
 114
 115 void 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
 129 void 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
 173 void 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
 180 void 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
 203 int 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
 211 class gen_vertices {
 212 public:
 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
 222 class gen_edges {
 223 public:
 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
 236 class reset_vertices {
 237 public:
 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
 248 void 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
 292 void ReleaseGraph() {
 293     delete []locks;
 294 }
 295
 296 void 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
 311 int 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 }