c_graph_ops_example.cpp

static inline int encode_distance(int start, int distance) { return ((start + 1) << 4) + distance; }

static inline int decode_distance(int x) { return x & ((1 << 4) - 1); }

static inline int node_visited(int x, int v) { return (x >> 4) == v + 1; }

template <typename Graph>
static void extend_closed_neighborhood(Graph const& G, Graph& H, int v, int radius, int buffer[]) {
  assert(0 <= v && v < G.number_of_nodes());
  assert(radius >= 0);

  // limitations
  if (G.number_of_nodes() > 100000000) abort();
  if (radius > 15) abort();

  std::queue<int> q;
  q.push(v);
  buffer[v] = encode_distance(v, 0);

  // run BFS
  while (!q.empty()) {
    auto parent = q.front();
    q.pop();

    for (auto u : G.neighbors(parent)) {
      if (!node_visited(buffer[u], v)) {
        // advance the frontier
        auto distance = decode_distance(buffer[parent]) + 1;
        buffer[u] = encode_distance(v, distance);

        // do not call `add_edge()` as it requires inter-thread locks
        H.edges_[v].push_back(u);
        if (distance < radius) q.push(u);
      }
    }
  }
}

template <typename Graph>
Graph distance_closure(Graph const& G, int radius, int num_threads) {
  assert(radius >= 0);
  assert(num_threads >= 1);
  typedef std::pair<int, int> PI;

  // special cases
  if (radius == 1) return G;

  auto n = G.number_of_nodes();
  Graph H(n);
  if (radius == 0) return H;

#pragma omp parallel num_threads(num_threads)
  {
    // create thread-local buffer
    int* buffer = (int*)std::calloc(n, sizeof(int));

    // main logic
#pragma omp for schedule(auto)
    for (int i = 0; i < n; ++i) { extend_closed_neighborhood(G, H, i, radius, buffer); }

    // free buffer
    std::free(buffer);
  }

  return H;
}

//
// Old version
//

// template <typename T>
// inline bool contains(std::unordered_set<T> const& s, T const& x) {
//   return s.find(x) != s.end();
// }

// std::unordered_set<int> get_closed_neighborhood(CIntGraph const& G, int v, int radius) {
//   assert(0 <= v && v < G.number_of_nodes());
//   assert(radius >= 0);
//   typedef std::pair<int, int> PI;

//   std::unordered_set<int> visited;
//   std::queue<PI> q;
//   q.push(std::make_pair(radius, v));

//   // run BFS
//   while (!q.empty()) {
//     auto p = q.front();
//     q.pop();
//     auto x = p.first;
//     auto y = p.second;

//     if (contains(visited, y)) continue;
//     visited.insert(y);
//     if (x == 0) continue;

//     for (auto u : G.neighbors(y)) {
//       if (!contains(visited, u)) {
//         // advance the frontier
//         q.push(std::make_pair(x - 1, u));
//       }
//     }
//   }
//   return visited;
// }

// CIntGraph distance_closure(CIntGraph const& G, int radius, int num_threads = 1) {
//   assert(radius >= 0);
//   typedef std::pair<int, int> PI;

//   // special cases
//   if (radius == 1) return G;

//   auto n = G.number_of_nodes();
//   CIntGraph H(n);
//   if (radius == 0) return H;

//   // main logic
//   if (num_threads <= 1) {
//     // single-threaded
//     for (int i = 0; i < n; ++i) {
//       for (auto j : get_closed_neighborhood(G, i, radius)) {
//         if (i < j) { H.add_edge(i, j); }
//       }
//     }
//   } else {
//     // multi-threaded
//     std::vector<PI>* results[num_threads];
//     for (int i = 0; i < num_threads; ++i) { results[i] = new std::vector<PI>; }

// #pragma omp parallel for schedule(auto) num_threads(num_threads)
//     for (int i = 0; i < n; ++i) {
//       for (auto j : get_closed_neighborhood(G, i, radius)) {
//         if (i < j) { results[omp_get_thread_num()]->push_back(std::make_pair(i, j)); }
//       }
//     }

//     // collect results
//     for (int i = 0; i < num_threads; ++i) {
//       for (auto& p : *results[i]) { H.add_edge(p.first, p.second); }
//       delete results[i];
//     }
//   }
//   return H;
// }