/* SSSP.cpp Single-source shortest path finder. Includes a (sequential) implementation of Dijkstra's algorithm, which is O((m + n) log n). Also includes a (sequential) implementation of Delta stepping. You need to create a parallel version of this. (c) Michael L. Scott, 2026. Based heavily on a Java version by the same author, written in 2022, translated to C++ with the help of ChatCPT, and subsequently modified and enhanced. For use by students in CSC 2/458 at the University of Rochester, during the Spring 2026 term. All other use requires written permission of the author. */ #include #include #include // for priority_queue #include #include #include #include #include #include // for getopt() using namespace std; // ============================================================ // Command-line configuration // ============================================================ struct Config { int n = 1000000; // default number of vertices double geom = 1.0; // default degree of geometric reality // 0 means random edge weight; // 1 means fully geometric distance int degree = 5; // expected number of neighbors per vertex // (near the middle of the graph) uint64_t seed = 0; // default random number seed unsigned int num_threads = 0; // zero means use Dijkstra's alg; // positive means use Delta stepping bool verbose = false; bool elog = false; bool wlog = false; }; void usage() { cerr << "usage:\n" "\t-n \n" "\t-d \n" "\t (for vertices near the middle of large graphs)\n" "\t-g \n" "\t (real number between 0 and 1)\n" "\t-s \n" "\t-t \n" "\t (0 means use Dijkstra's algorithm on one thread)\n" "\t-v print verbose diagnostics\n" "\t-e print log of edges (de)selected\n" "\t-w print end-of-run vertex path weights\n" "\t-h print this message\n"; exit(0); } Config parse_args(int argc, char** argv) { Config cfg; int ch; while ((ch = getopt(argc, argv, "n:d:g:s:t:vewh")) != -1) { try { switch(static_cast(ch)) { case 'n': cfg.n = stoi(optarg); break; case 'd': cfg.degree = stoi(optarg); break; case 'g': cfg.geom = stod(optarg); break; case 's': cfg.seed = stoull(optarg); break; case 't': cfg.num_threads = stoul(optarg); break; case 'v': cfg.verbose = true; break; case 'e': cfg.elog = true; break; case 'w': cfg.wlog = true; break; case 'h': // fall through case '?': usage(); } } catch (std::exception const& ex) { cerr << "invalid argument to -" << static_cast(ch) << ": " << optarg << endl; usage(); } } return cfg; } // ============================================================ // Graph representation // ============================================================ struct Edge; constexpr auto UNREACHABLE = numeric_limits::max(); struct Vertex { int x; int y; vector neighbors; uint64_t dist = UNREACHABLE; Edge* pred = nullptr; Vertex(int _x, int _y) : x(_x), y(_y) {}; }; struct VertexHash { size_t operator()(const Vertex* v) const noexcept { return static_cast(v->x) ^ (static_cast(v->y) << 1); } }; struct Edge { Vertex* a; // vertices are in arbitrary order Vertex* b; int weight; Vertex* other(Vertex* v) const { return v == a ? b : a; } Edge(Vertex* _a, Vertex* _b, int w) : a(_a), b(_b), weight(w) {}; }; class Graph { public: static constexpr int MAX_COORD = 1024 * 1024 * 256; Graph(Config cfg) : cfg_(cfg), prng_(cfg.seed) { cout << cfg.n << " vertices, "; if (cfg.num_threads) { cout << "Delta (" << cfg.num_threads << ")"; } else { cout << "Dijkstra"; } cout << ", degree " << cfg.degree << ", geometry " << cfg.geom << ", seed " << cfg.seed << endl; if (cfg_.verbose) { cout << "initializing vertices... " << flush; } vertices_.reserve(cfg_.n); uniform_int_distribution coord_dist(0, MAX_COORD - 1); unordered_set vertex_hash; // Create vertices with unique coordinates while (static_cast(vertices_.size()) < cfg_.n) { auto v = make_unique(coord_dist(prng_), coord_dist(prng_)); // v will be a unique_ptr, with ownership of the new Vertex if (vertex_hash.insert(v.get()).second) { // not a duplicate vertices_.push_back(std::move(v)); // transfers ownership of Vertex from v to vertices_ } } vertices_[0]->dist = 0; // source vertex buildEdges(); if (cfg_.verbose) { cout << endl << "done"; if (cfg_.num_threads) { // get this message out before timing begins cout << "; delta = " << MAX_COORD / cfg_.degree; } cout << endl << flush; } } // ------------------------------------------------------------ // Dijkstra's algorithm // ------------------------------------------------------------ // Compare vertices by distance from root. Use greater-than so pq sinks // heavy nodes, leaving the minimum at the root. struct ShorterDistance { constexpr bool operator()(const Vertex* lhs, const Vertex* rhs) const { return lhs->dist > rhs->dist; } }; void dijkstraSolve() { priority_queue, ShorterDistance> pq; Vertex* v = vertices_[0].get(); // returns pointer to vertex // but leaves ownership with the vertices_ vector while (true) { for (auto* e : v->neighbors) { auto* o = e->other(v); uint64_t alt = v->dist + e->weight; if (alt < o->dist) { o->dist = alt; o->pred = e; pq.push(o); } } if (pq.empty()) break; v = pq.top(); pq.pop(); selectEdge(v->pred); } } // ------------------------------------------------------------ // Delta stepping (currently sequential) // ------------------------------------------------------------ void deltaSolve() { const int NUM_BUCKETS = 2 * cfg_.degree; const int delta = MAX_COORD / cfg_.degree; // All buckets, together, cover a range of 2 * MAX_COORD, // which is larger than the weight of any edge, so a relaxation // will never wrap all the way around the array. vector> buckets(NUM_BUCKETS); // each bucket is a hash table, with (expected) const-time // insert and remove (erase) buckets[0].insert(vertices_[0].get()); // local function (lambda) // reconsider path to o: is it easier to reach via e? auto relax = [&](Vertex* o, Edge* e) { Vertex* v = e->other(o); uint64_t alt = v->dist + e->weight; if (alt < o->dist) { int oldb = (o->dist / delta) % NUM_BUCKETS; int newb = (alt / delta) % NUM_BUCKETS; buckets[oldb].erase(o); o->dist = alt; if (o->pred) unselectEdge(o->pred); o->pred = e; selectEdge(e); buckets[newb].insert(o); } }; int i = 0; vector> heavy_edges; // for vertices guaranteed to relax into a later bucket while (true) { while (!buckets[i].empty()) { unordered_set contents = std::move(buckets[i]); buckets[i].clear(); for (Vertex* v : contents) { for (Edge* e : v->neighbors) { if (e->weight <= delta) { // lightweight relaxation; // might end up in buckets[i] again relax(e->other(v), e); } else { heavy_edges.push_back({e->other(v), e}); } } } } // Nothing will ever go back in bucket i at this point. // Now consider remembered edges of weight > delta: for (auto [o, e] : heavy_edges) { relax(o, e); } heavy_edges.clear(); // Find next nonempty bucket: int j = i; do { j = (j + 1) % NUM_BUCKETS; } while (j != i && buckets[j].empty()); if (j == i) { // Cycled all the way around; we're done break; // while (true) lop } i = j; } } void report() { if (cfg_.wlog) { cout << "final vertices (x y dist):" << endl; for (auto& v : vertices_) { cout << "\t" << setw(10) << v->x << setw(10) << v->y << setw(22) << v->dist << endl; } } if (cfg_.verbose && cfg_.num_threads) { cout << deselections << " deselections" << endl; } } private: // ------------------------------------------------------------ // Checkerboard-based edge construction // ------------------------------------------------------------ void buildEdges() { // As a heuristic, I want to connect each vertex to about 1/4 of // its geometrically nearby vertices. So I want to choose // neighbors from a region containing about 4*degree vertices. // I divide the plane into a k x k grid, such that a 3x3 subset // has about the right number of vertices from which to choose. const int k = static_cast( sqrt(static_cast(cfg_.n) / cfg_.degree) * 3.0 / 2.0); const int sw = static_cast( // width of checkerboard square ceil(static_cast(MAX_COORD) / k)); if (cfg_.verbose) { cout << "and edges (" << k << " sq grid, " << setprecision(1) << fixed << (double) cfg_.n / k / k << " pts per cell) ... " << flush; } vector>> board( k, vector>(k)); // Place vertices into checkerboard cells for (auto& v : vertices_) { int i = v->x / sw; int j = v->y / sw; board[i][j].push_back(v.get()); } uniform_int_distribution coin(0, 3); uniform_int_distribution rand_weight(0, MAX_COORD * 2); for (auto& v : vertices_) { int xb = v->x / sw; int yb = v->y / sw; int xl, xh, yl, yh; if (k < 3) { xl = yl = 0; xh = yh = k - 1; } else { xl = (xb == 0) ? 0 : (xb == k - 1 ? k - 3 : xb - 1); xh = (xb == 0) ? 2 : (xb == k - 1 ? k - 1 : xb + 1); yl = (yb == 0) ? 0 : (yb == k - 1 ? k - 3 : yb - 1); yh = (yb == 0) ? 2 : (yb == k - 1 ? k - 1 : yb + 1); } VertexHash h; // for deterministic choice among vertices across program runs for (int i = xl; i <= xh; ++i) { for (int j = yl; j <= yh; ++j) { for (Vertex* u : board[i][j]) { if (h(v.get()) < h(u) // only choose edge from one end -- // avoid self-loops and doubled edges && coin(prng_) == 0) { // 1/4 chance // invent a weight: double dx = v->x - u->x; double dy = v->y - u->y; int dist = static_cast(sqrt(dx * dx + dy * dy)); int w = static_cast( cfg_.geom * dist + (1.0 - cfg_.geom) * rand_weight(prng_)); // create an edge between u and v auto e = make_unique(v.get(), u, w); v->neighbors.push_back(e.get()); u->neighbors.push_back(e.get()); edges_.push_back(std::move(e)); } } } } } } void selectEdge(Edge* e) { if (cfg_.elog) { cout << "selected " << setw(10) << e->a->x << setw(10) << e->a->y << " " << setw(10) << e->b->x << setw(10) << e->b->y << " " << setw(10) << e->weight << endl; } } void unselectEdge(Edge* e) { if (cfg_.elog) { cout << "unselected " << setw(10) << e->b->x << setw(10) << e->b->y << " " << setw(10) << e->weight << endl; } deselections++; } Config cfg_; mt19937_64 prng_; // pseudo-random number generator int deselections = 0; vector> vertices_; // owns vertices vector> edges_; // owns edges }; // ============================================================ // main // ============================================================ int main(int argc, char** argv) { auto cfg = parse_args(argc, argv); Graph g(cfg); // Don't start the clock until the graph has been constructed. auto start = chrono::steady_clock::now(); if (cfg.num_threads == 0) { g.dijkstraSolve(); } else { g.deltaSolve(); } auto end = chrono::steady_clock::now(); chrono::duration elapsed = end - start; g.report(); // provide any desired end-of-run information // elapsed time is meaningless if we've been logging edges if (!cfg.elog) { cout << "elapsed time: " << elapsed.count() << " seconds" << endl; } }