oneAPIを使ったFPGA上でのグラフ探索アルゴリズムの実装

#1 はじめに
この記事は、インテル® FPGA Advent Calendar 2021の23日目の記事です。
oneAPIとはインテル®社が提供しているクロスアーキテクチャプログラミングモデルでありC++ベースの言語(DPC++)を用いて複数のアーキテクチャ向けソフトウェアを容易に開発できるものになってます。
今回はそのDPC++を使って、FPGAで動作する簡単なグラフ上のアルゴリズムである幅優先探索を実装しました。幅優先探索（Breadth-First Search, BFS）とは、グラフ上の始点殻から全てのノードまでの最短経路を見つけることができ、多くの応用で用いられています。また、スーパーコンピュータを対象としたベンチマークGraph 500 でも採用されているアルゴリズムです。

#2 動作環境

• OS : CentOS Linux release 7.9.2009 (Core)
• CPU : Intel® Xeon® Gold 6242 CPU @ 2.80GHz
• FPGA : Intel® Stratix® 10 GX FPGA

#3 グラフについて
グラフとはデータの世界では下の図のようなものを言います

0から8がノードと呼ばれるものです。
図からわかりますが線でつながっているノード同士があります。これはそのノード同士が隣接していることを表しています。
この線のことをエッジと言います。
(今回は簡略化のために扱うグラフは無向グラフにします)
このようなグラフは一般的に以下のようなtxtファイルやcsvファイルとして保存されています。

data.txt

0 1
0 4 
0 2
・　　・
・　　・

data.txtの1行目はノード0とノード1が隣接しているという意味で2行目はノード0とノード4が隣接してるという意味になってます。

グラフの表現

グラフの表現方法には様々な種類がありますが今回は隣接行列を用いることにします。

隣接行列

隣接行列とは各ノードがどのノードに対してエッジを持っているかを表している行列だと思ってください。
例えば上のグラフではノード0はノード1,2,4とつながっているので行列の0行目の1,2,4列目に1がそれ以外には0が格納されます。(下の行列のように)
つまり、隣接行列はノード数×ノード数の行列になります。

\begin{pmatrix}
0 & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 1 \\
1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
\vdots  & & & & \vdots & & & & \vdots
\end{pmatrix}

#4 幅優先探索
ここからは幅優先探索(以下BFS)の具体的なコードについて説明します。
##BFSのアルゴリズム
BFSの具体的な処理について説明します。
まずは探索を始める根ノードを決めます。
今回はノード0を根ノードにするのでまずノード0が訪問済みとなります。
次に、ノード0の隣接ノード1,2,4が訪問済みとなり、その次はノード1,2,4の隣接ノードが訪問済みとなり...と隣接ノードをたどっていくことになります。BFSの処理が終了するのはすべてのノードが訪問済みになったときです。
図を使って具体的に説明します。

まず、訪問するノードを格納する配列(この配列のことをフロンティアテーブルと言います)にノード0が格納される。

次にノード0の隣接ノードであるノード1,2,4がフロンティアテーブルに格納されます。
このとき、ノード0からの距離を格納する配列(dist)にノード0からノード1,2,4までの距離1が格納されます。
この一連の動作を1ステップと言います。なのでdistに格納する距離はステップ数と同じになります。

この動作をすべてのノードが訪問済みとなるまで繰り返します。

#実装
これから出てくる定数、配列の意味は次の通りです。

int edge_array[node_count][node_count] = {};  // 隣接リスト
const size_t node_count = 9;         // ノード数
int now_frontier[node_count] = {};   // 現在のレベルのフロンティア
int next_frontier[node_count] = {};  // 次のレベルのフロンティア
int dist[node_count] = {};           // ルートからの距離
int step = 1;                       　　// ステップ数
int converged[0] = {};               // 探索の終了条件

edge_arrayは対象のグラフの隣接行列が格納されています。
now_frontierは現在のステップでのフロンティア、next_frontierは次のステップでのフロンティアとなります。
distには根ノードからの各ノードへの最短距離が格納されることになります。

##①BFSを実装するソースコード

第一段階(①)

    q.submit([&](handler &h) {
      accessor edge_array(edge_buf, h, read_only);
      accessor now_frontier(now_frontier_buf, h, read_write);
      accessor next_frontier(next_frontier_buf, h, read_write);

      h.single_task([=]() {
        for (int i = 0; i < node_count; i++) {
          if (now_frontier[i] == 1) {
            for (size_t j = 0; j < node_count; j++) {
              if (edge_array[i][j] == 1) {
                next_frontier[j] = 1;
              };
            };
          };
        };
      });
    });
    q.wait();

実際のBFSは６行目のfor文からです。現在のステップでのフロンティアに格納されているノードに対して隣接しているノードを見つけます。(9行目のif文)
ここで見つかった隣接ノードは次のステップでフロンティアとして使用します。

##②距離計算
距離計算をするソースコードです。

第二段階(②)

    q.submit([&](handler &h) {
      accessor now_frontier(now_frontier_buf, h, read_only);
      accessor next_frontier(next_frontier_buf, h, read_only);
      accessor dist(dist_buf, h, write_only);
      accessor converged(converged_buf, h, write_only);

      h.single_task([=]() {
        for (int i = 0; i < node_count; i++) {
          if (now_frontier[i] == 0 && next_frontier[i] == 1) {
            dist[i] = level;
          };
        };
      });
    });
    q.wait();

距離を格納する条件は現在のステップでそのノードが未訪問であり、かつ次のステップでそのノードが訪問予定で有ることです。
言い換えるとnow_frontierでは0がnext_frontierでは1が格納されていることが条件となります。
この条件を書いているのが7行目のif文になります。

##③最後に
next_frontierをnow_frontierにする必要があるのでこの2つの配列要素を交換します。
これはステップnのnext_frontierはステップn+1でのnow_frontierとなるからです。
また、next_frontierとnow_frontierが入れ替え終わったあとにステップ数を+１します。

第三段階(③)

    q.submit([&](handler &h) {
      accessor now_frontier(now_frontier_buf, h, read_write);
      accessor next_frontier(next_frontier_buf, h, read_only);

      h.single_task([=]() {
        for (int i = 0; i < node_count; i++) {
          now_frontier[i] = next_frontier[i];  //フロンティアをswap
        };
      });
    });
    q.wait();
    step++;

#全体のソースコードの流れ
先程述べた①から③の処理をすべてのノードが訪問済みになるまで繰り返します。
これは言い換えると、距離計算が行われなかったときに処理を終了するということになります。
なので①から③のソースコードをwhile文でくくり②のソースコードを以下のように変更します。

    q.submit([&](handler &h) {
      accessor now_frontier(now_frontier_buf, h, read_only);
      accessor next_frontier(next_frontier_buf, h, read_only);
      accessor dist(dist_buf, h, write_only);
      accessor converged(converged_buf, h, write_only);

      h.single_task([=]() {
        for (int i = 0; i < node_count; i++) {
          if (now_frontier[i] == 0 && next_frontier[i] == 1) {
            dist[i] = level;
            converged[0] = 0;  //距離が生成されたら探索を続ける
          };
        };
      });
    });
    q.wait();

全体の流れ

while(!converged[0]) {
　　 converged[0] = 1;
  ①
　　　　②
　　　　③
}

しかし、この構成ではBFSは正常に動作しません。
それはwhile文を回す条件のcenverged配列で問題があるからです。
converged配列が②で更新されたときこの更新はFPGA上で更新されているだけであり、CPU上では反映されていない可能性があるからです。反映されていないとwhile文を回す条件を満たしていないことになるのですぐにwhile文から抜け出してしまい、正常にBFSが動作しない事になってしまいます。
なので、converged配列を変更したあとにconverged配列をCPU、FPGAで同期させる必要があります。
そのコードが以下の2つです。

④(CPUからFPGAへの同期)

    q.submit([&](handler &h) {
      accessor converged_sync(converged_buf, h, write_only);
      h.copy(converged, converged_sync);
    });

⑤(FPGAからCPUへの同期)

    q.submit([&](handler &h) {
      accessor converged_sync(converged_buf, h, read_only);
      h.copy(converged_sync, converged);
    });

この④、⑤をコードに追加することでBFSが正常に起動するようになります。

全体の流れ(修正版)

while(!converged[0]) {
　　 converged[0] = 1;
  ①
　　　　②
　　　　④
　　　　③
　　　　⑤
}

これでBFSを実装したプログラムの解説を終わりにします。

#実行結果
最後に実行結果を貼っておきます。

Running on device: pac_s10 : Intel PAC Platform (pac_f000000)
Node size: 9
dist[0]:  = 0
dist[1]:  = 2
dist[2]:  = 2
...
dist[8]:  = 3
bfs successfully completed on device.