はじめに

MDACPのGPGPUへの対応を行い一応完了したので、作業の過程とベンチマークの結果をここにまとめておく。ソースコードはここにある。

MDACPとは

物性研究所の渡辺宙志助教らによって開発された並列分子動力学シミュレーションコード。version 1ではflat MPI、version 2ではOpenMPとMPIのハイブリッドでの実行ができるようになっている。本家のホームページはここになる。今回はversion 2のGPGPU化を行なった。

MDACPで採用されている並列化

基本的には空間を均一にメッシュで切って領域を分割したのち、分割した小領域にスレッドを対応させている。
ここら辺についてはslideshareに情報が公開されていたので、こちらを参考にした。

移植する上で割と重要だった情報は
1. mdunitがスレッドに対応している。mdunitが分割された領域の最小単位に相当する。
2. mdmanagerがプロセスに対応している。
3. variablesクラス内にあるGetTotalParticleNumberとGetParticleNumberの違い。

1と2についてはコードをみてすぐにわかったのだが、3についてはslideshareの情報がなければ正直わからなかった。

GPGPU化の方針

まず、どの部分をオフロードするかを決める。分子動力学シミュレーションの手順は基本的には

  1. 初期配置と速度の生成 (これは別プログラムとして独立させることが多い。)
  2. Verlet listの構築
  3. 力の計算
  4. 位置と運動量のアップデート
  5. Verlet listが有効であるかどうかを確認 (無効ならば2に飛ぶ)
  6. 目的ステップ数未満であれば、3に飛ぶ
  7. 終了

のようになる。ここでホットスポットは3の力の計算になる。よってこの部分をGPUで計算させることにする。
GPUを使った分子動力学シミュレーションコードの中には、GPUで計算させている間、CPUが遊んでしまう実装と、GPUで計算させている間にCPUで別の処理をさせる実装(例えば、GROMACSなど)の二種類がある。今回はできるだけCPUを遊ばせないように、後者の実装を採用した。

GPGPU力計算用のカーネルとしては自作のこのコードを用いた。Pascal向け最速コードを採用した。理由としてはKepler向け最速コードがVerlet Listのデータレイアウトを元のmdacpで採用されているデータレイアウトから変更する必要があるから。

GPGPUによる力計算の流れとしては
1. Verletリストを構築したらそのデータを必要な分だけGPUへ転送する。
2. 力計算を行う前に、粒子の座標データと運動量データをGPUで仕事させる分だけCPU -> GPUへ転送。
3. GPUで力を計算させる。
4. 力計算が終わった後で運動量データだけGPU -> CPUへ転送
5. 2 ~ 4を実行させながら、CPU側でも力の計算を実行させる。
6. CPUとGPUで同期をとる。

CPUとGPUで担当させる仕事の割合をどう分散させるかについて。実行時にオートチューニングするようなことも考えたが、相当実装が大変そうだったので、ひとまず諦めた。

また、一つのGPUに対して一つのプロセスを割り当てることにする。(つまり、一つのmdmanagerが一つのGPUの計算を管理する。)結果として、複数のmdunitの計算を一つのGPUで行わせる。

作業内容

作業内容について、一部抜粋して以下に記す。一部GPGPUとは関係ない部分も含んでいる

cmakeの導入

もともとmakefile.optというものにコンパイラやコンパイルオプションを書き込む形式になっていたが、これをcmakeに変更した。cmakeの簡単な使い方や、CUDAをコンパイルするときの設定については、ごく簡単なcmakeの使い方CUDAもcmakeでコンパイルするを参考のこと。

二つの記事などを参考にしながら、CMakeLists.txtを作成したのち、src include ディレクトリに必要なファイルを放り込んでコンパイルできるようにした。

GPUとCPUのデータ転送を行うクラス (CudaPtrとCudaPtr2D)を追加

cuda_ptr.hcuda_ptr2d.hになる。
CudaPtrはデバイスとホストのデータを保持するためのクラスで、CudaPtr2Dは用途は同じだが座標と運動量データを保持するための専用クラスになっている。

二つクラスを分ける必要性ははっきり言ってない。
なぜこんなことをしたのかというと、mdacpのコードでは座標と運動量データが

class Variables {
public:
...
  __attribute__((aligned(64))) double q[N][D];
  __attribute__((aligned(64))) double p[N][D];
...
};

のように、二次元配列として宣言されていたから。
CUDAのカーネル関数で二次元配列を使うことはできなくはない。だが、少し面倒なので、デバイスでは一次元配列としてデータを扱いたかった。そのため、ホスト側で二次元配列のデータをデバイス側で一次元配列として扱うCudaPtr2Dを専用で用意した。

二つのクラスにはホストデバイス間の通信を行うメンバー(ブロッキングとノンブロッキングの双方を用意)を用意している。例えば、

  void Host2Dev(const int beg,
                const int count) {
    checkCudaErrors(cudaMemcpy(dev_ptr_ + beg,
                               host_ptr_ + beg,
                               count * sizeof(T),
                               cudaMemcpyHostToDevice));
  }
  void Host2Dev(void) {this->Host2Dev(0, size_);}
  void Host2DevAsync(const int beg,
                     const int count,
                     cudaStream_t strm = 0) {
    checkCudaErrors(cudaMemcpyAsync(dev_ptr_ + beg,
                                    host_ptr_ + beg,
                                    count * sizeof(T),
                                    cudaMemcpyHostToDevice,
                                    strm));
  }

である。Asyncなしはブロッキング版でAsyncありはノンブロッキング版である。

時間計測クラスにGPUのkernel関数の実行時間を測定するメソッドを追加。

元のmdacpのコードにStopWatchクラスが用意されていたが、GPUのカーネルの実行時間を取るためにこのStopWatchクラスは使えない。
cudaDeviceSynchronizeした上での時間計測するのであれば良いのだけども、それだとGPUの処理が終わるまでCPUの処理がストップしてしまうので、CPUとGPUで同時に力計算を行うといったことができなくなってしまう。

そういうときにはcudaEventRecordを使って、次のように時間計測を行えば良い。

// GPU execution
float elapsed_time = 0.0;

cudaEventRecord(start);
kernel<<<gr_size, tb_size>>>(dev_ptr);
cudaEventRecord(stop);

// CPU execution
do_something_on_cpu();

// CPU GPU synchronization
cudaDeviceSynchronize();

cudaEventSynchronize(stop);
float elapsed_time = 0.0;
cudaEventElapsedTime(&elapsed_time, start, stop);

ここで、cudaEventRecordを呼んだ時点で同期処理は入らず、do_something_on_cpu()の実行される。
elapsed_timeにはkernel関数の実行時間が格納される。

cudaEvent_tをStopWatchに実装し、GPUでのkernelの実行時間を取れるように変更した。使い方は元のStopWatchクラスと同じ。

  // gpu force calc
  swForce_gpu.Start();
  for (int i = 0; i < num_threads; i++) { // GPU
    pn_gpu[i] = int(mdv[i]->GetParticleNumber() * WORK_BALANCE);
    mdv[i]->CalculateForceGPU(pn_gpu[i], strms[i]);
  }
  swForce_gpu.Stop();

  // cpu force calc
  swForce_cpu.Start();
  #pragma omp parallel for schedule(static)
  for (int i = 0; i < num_threads; i++) { // CPU
    mdv[i]->CalculateForceCPU(pn_gpu[i]);
  }
  swForce_cpu.Stop();

  // sync
  checkCudaErrors(cudaDeviceSynchronize());
  swForce_gpu.Record();

実装についてはstopwatch.hを参照。

GPUデバイスの情報を表示するための関数(device_query_all)を追加。

ノード内に刺さっている全てのGPUの情報を表示するための関数。cuda runtime APIを順次呼ぶだけのもので、特に工夫はない。device_info.cc

cudaStreamを使ったCPU <-> GPUの転送時間の隠蔽

1GPUに複数のmdunitを割り当てるが、各mdunitの力計算は独立に行えるので、全てデフォルトストリームで計算させるのは勿体無い。
そこで、mdunitごとにcudaStreamを作成し、非同期に計算させることにする。
こうすることの最大のメリットはCPU <-> GPUの転送時間の隠蔽ができることである。
nvvpのタイムラインを見た方がわかりやすいと思うが、

Screenshot from 2017-04-20 13:53:12.png

のようになる。茶色の部分がCPU <-> GPUの転送時間に対応し、緑色の部分が力計算の部分に対応する。CPU <-> GPUの転送と同時にGPUでの計算が行われているので、CPU <-> GPUの転送時間がほぼ見えなくなっている。mdunit8つに対してcudaStreamを8つ作っており、CPU <-> GPUの転送時間は本来の転送時間の1/8になっていることがわかる。

ベンチマーク

ベンチマークは物性研究所システムBを用いて行った。
実行環境は以下の通り

CUDA version 7.0
gcc version 4.8.5
SGI MPT version 2.14

また、インプットファイルは次の通り

inpug_gpu.cfg
Mode=Benchmark
Density=0.712
ThermalizeLoop=150
UnitLength=45
TimeStep=0.01
TotalLoop=1000
InitialVelocity=0.9

物性研究所システムBを使って64GPUs、512CPUcoresで計算させた時の結果になる。
粒子数は33219072個でおよそ3千万粒子。1000ステップ計算するのにかかった時間を表している。

  • GPUありの場合 33.93638 [SEC]
  • GPUなしの場合 61.47893 [SEC]

実行時間としては2倍弱の高速化に止まっている。

計算時間の内訳

計算時間の内訳は次のようになる。

comm [sec] force [sec] pair [sec] other [sec]
CPU + GPU 17.6621 11.3571 2.66857 2.24861
CPU 18.1391 39.5660 2.68539 1.08844

もともとMPI通信時間が大きいみたいで、これをもう少し時間削減して、力計算部分の時間の割合を増やすことができれば、GPGPU化の効果がもう少し出てくるはずである。

まとめ

mdacpのGPGPU移植を行い、元の実装から2倍弱の高速化に成功した。
擬似flat MPI実装からのGPGPU移植は正直無理なんじゃないかと思っていたが、やってみると意外とできるもんだなぁ(小並感)。

今後

AVX512対応?

編集履歴

  • 表のデータに誤りがあったので修正。