GAP Benchmarkの各プログラムの構成について (sssp.ccのソースコード)

GAP Benchmarkの各プログラムの中身を調べている。sssp.ccについて調べた。

sssp: Single source shortest path

SSSPのコードをもうちょっと読み解いていこうと思う。

github.com

重要な変数：

• distベクトル：各ノードまでの最短距離を保持する。つまり、始点からノードiまでの最短距離をdist[i]として表現する。 最初はkDistInfという定数でほぼ無限大を表現している。その後、アルゴリズムが進行するとRelaxEdge()と言って、各ノードまでの最短距離が更新される。
• frontierベクトル：協会ノードを保持する。現在のイタレーションで訪問すべきノードの集合を表現する。
• local_binsベクトル：エッジのRelaxationの結果として得られた新しい最短距離に基づいて、ノードと適切のバケットに配置する。
  • 具体的には、新しい最短距離new_distが発見されると、対応するノードはnew_distに基づいて適切なバケット（local_bins[new_dist/delta]）に配置される。
  • そして、アルゴリズムの各イテレーションでは現在のバケット（local_bins[curr_bin_index]）に含まれるノードのエッジがRelaxationされる。
  • このコードの実装では各スレッドが自身のlocal_binsを持つようになっている。
    • これにより、異なるスレッド間でlocal_binsへのアクセスが競合することを防ぎ、並列処理の効率を高める。
  • 各イタレーションの終わりには、local_binsからfrontierへとノードがコピーされ、次のfrontier（つまり次のイタレーションで訪問すべきノードセット）が形成される。

pvector<WeightT> DeltaStep(const WGraph &g, NodeID source, WeightT delta) {
  Timer t;
  pvector<WeightT> dist(g.num_nodes(), kDistInf);
  dist[source] = 0;
  pvector<NodeID> frontier(g.num_edges_directed());
  // two element arrays for double buffering curr=iter&1, next=(iter+1)&1
  size_t shared_indexes[2] = {0, kMaxBin};
  size_t frontier_tails[2] = {1, 0};
  frontier[0] = source;
  t.Start();
  #pragma omp parallel
  {
    vector<vector<NodeID> > local_bins(0);
    size_t iter = 0;
    while (shared_indexes[iter&1] != kMaxBin) {
      size_t &curr_bin_index = shared_indexes[iter&1];
      size_t &next_bin_index = shared_indexes[(iter+1)&1];
      size_t &curr_frontier_tail = frontier_tails[iter&1];
      size_t &next_frontier_tail = frontier_tails[(iter+1)&1];
      #pragma omp for nowait schedule(dynamic, 64)
      for (size_t i=0; i < curr_frontier_tail; i++) {
        NodeID u = frontier[i];
        if (dist[u] >= delta * static_cast<WeightT>(curr_bin_index))
          RelaxEdges(g, u, delta, dist, local_bins);
      }
      while (curr_bin_index < local_bins.size() &&
             !local_bins[curr_bin_index].empty() &&
             local_bins[curr_bin_index].size() < kBinSizeThreshold) {
        vector<NodeID> curr_bin_copy = local_bins[curr_bin_index];
        local_bins[curr_bin_index].resize(0);
        for (NodeID u : curr_bin_copy)
          RelaxEdges(g, u, delta, dist, local_bins);
      }
      for (size_t i=curr_bin_index; i < local_bins.size(); i++) {
        if (!local_bins[i].empty()) {
          #pragma omp critical
          next_bin_index = min(next_bin_index, i);
          break;
        }
      }
      #pragma omp barrier
      #pragma omp single nowait
      {
        t.Stop();
        PrintStep(curr_bin_index, t.Millisecs(), curr_frontier_tail);
        t.Start();
        curr_bin_index = kMaxBin;
        curr_frontier_tail = 0;
      }
      if (next_bin_index < local_bins.size()) {
        size_t copy_start = fetch_and_add(next_frontier_tail,
                                          local_bins[next_bin_index].size());
        copy(local_bins[next_bin_index].begin(),
             local_bins[next_bin_index].end(), frontier.data() + copy_start);
        local_bins[next_bin_index].resize(0);
      }
      iter++;
      #pragma omp barrier
    }
    #pragma omp single
    cout << "took " << iter << " iterations" << endl;
  }
  return dist;
}