LJ力計算のGPGPU化 作用反作用の法則を使わないやり方がベストなのか?

はじめに

(2016/1/27追記)
致命的な間違いをおかしており、この記事の結果が変わりました。結局チューニングするとギリギリで作用反作用使わないほうが高速化されました。後日修正した記事を書く予定です。

これまでLJ力計算のGPGPU化では作用反作用を無視して、atomicを使わないでやるのが最速と自分は思ってきました。過去の文献を調べてみても(かなり前の論文ですが)この論文では

If such features (atomicのこと) did exist, the increment of Fk (j粒子の力) inside the inner
loop still requires a scattered memory access pattern which would slow performance significantly.

と、「j粒子への力の書き戻しは書き込みパターンが不連続になるので、仮にatomicあったとしても性能出ないと思うよ。」と言っています。(この論文出たころはまだ浮動小数点のatomicAddはありませんでした。)

そのあとに出てきた散逸粒子動力学法のGPGPU化の論文を見てみると、一生懸命作用反作用を使わないで計算する方法が提案されています。(散逸粒子動力学の説明はしないですが、ランダム力を粒子対について作用反作用を満たすように与えないといけないので、作用反作用なしで計算しようとすると非常に複雑な手続きが必要になります。そのため、GPGPU化、MPI化が難しいとされてきました。)

青木先生らの講義資料でも使用は最少限にとどめるべきとあります。

しかし、GPUは変化し続けています。現在のGPUでもLJ力計算を行う際に作用反作用を使わないのが本当にベストなのかどうかは正直不明です。そこで、作用反作用を使う場合とそうでない場合とで性能がどの程度違ってくるのかを調べてみました。

以下「実行時間」と書いたときは
1. GPU<->CPUの転送の時間込み。
2. 関数を100回続けて実行させて、実行時間の総和をとったもの。
になります。

ソースコードはここにあります。

LJ力計算のGPGPU化

簡単にCPUのコードを例にとってLJ力計算の処理を説明します。

ref.cpp
void force_cpu_ref(void){
  const int pn =particle_number;
  for (int i=0; i<pn; i++) { // すべての粒子について
    const double qx_key = q[i][X];
    const double qy_key = q[i][Y];
    const double qz_key = q[i][Z];
    const int np = number_of_partners[i]; // i粒子の周りにいる粒子の数
    const int kp = pointer[i];
    for (int k=0; k<np; k++) { // i粒子の周りにいるすべてのj粒子について
      const int j = sorted_list[kp + k];
      double dx = q[j][X] - qx_key;
      double dy = q[j][Y] - qy_key;
      double dz = q[j][Z] - qz_key;
      double r2 = (dx*dx + dy*dy + dz*dz); // 距離を計算し
      if (r2 > CL2) continue; // cutoff 距離以上であれば以降の処理をスキップ
      double r6 = r2*r2*r2;
      double df = ((24.0*r6-48.0)/(r6*r6*r2))*dt;
      p[i][X] += df*dx; // 作用反作用の法則を使い、計算した力をたしこむ。
      p[i][Y] += df*dy;
      p[i][Z] += df*dz;
      p[j][X] -= df*dx;
      p[j][Y] -= df*dy;
      p[j][Z] -= df*dz;
    }
  }
}

二重ループになっていて、並列に処理できる箇所が二つあります。

この関数のGPGPU化でよくやられるのはi粒子の数だけGPUのスレッドを起動させて並列処理するやりかただと思います。ナイーブな実装ですがatomicAddを用いて

作用反作用使う
const auto tid = threadIdx.x + blockIdx.x * blockDim.x; // スレッドがi粒子一つに対応する。
if (tid < particle_number) {
  const auto qi = q[tid];
  const auto np = number_of_partners[tid];
  const auto kp = pointer[tid];

  for (int32_t k = 0; k < np; k++) { // i粒子の周りにいるj粒子についてのループ
    const auto j = sorted_list[kp + k];
    const auto dx = q[j].x - qi.x;
    const auto dy = q[j].y - qi.y;
    const auto dz = q[j].z - qi.z;
    const auto r2 = dx * dx + dy * dy + dz * dz;
    if (r2 > CL2) continue;
    const auto r6 = r2 * r2 * r2;
    const auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
    atomicAdd(&p[tid].x, df * dx);
    atomicAdd(&p[tid].y, df * dy);
    atomicAdd(&p[tid].z, df * dz);
    atomicAdd(&p[j].x, -df * dx);
    atomicAdd(&p[j].y, -df * dy);
    atomicAdd(&p[j].z, -df * dz);
  }
}

となります。

作用反作用を使わないやり方の場合には、j粒子の運動量の書き戻しが発生しないため、atomicAddがなくなり、下記のようになります。

作用反作用使わない
const auto tid = threadIdx.x + blockIdx.x * blockDim.x; // i粒子一つに対応する。
if (tid < particle_number) {
  const auto qi = q[tid];
  const auto np = number_of_partners[tid];
  const auto kp = pointer[tid];

  for (int32_t k = 0; k < np; k++) { // i粒子の周りにいるj粒子についてのループ
    const auto j = sorted_list[kp + k];
    const auto dx = q[j].x - qi.x;
    const auto dy = q[j].y - qi.y;
    const auto dz = q[j].z - qi.z;
    const auto r2 = dx * dx + dy * dy + dz * dz;
    if (r2 > CL2) continue;
    const auto r6 = r2 * r2 * r2;
    const auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
    p[tid].x += df * dx;
    p[tid].y += df * dy;
    p[tid].z += df * dz;
  }                                                                           
}

単精度の場合と倍精度の場合で上記二つの関数の実行時間を比較してみると以下のようになります。
GPU: Tesla K40t
粒子数:10万
密度: 1

実行時間 [s]
作用反作用あり (単精度) 0.885239
作用反作用なし (単精度) 1.361835
作用反作用あり (倍精度) 2.450472
作用反作用なし (倍精度) 1.996499

密度が0.5で粒子数が62500の時

実行時間 [s]
作用反作用あり (単精度) 0.224812
作用反作用なし (単精度) 0.379074
作用反作用あり (倍精度) 0.696441
作用反作用なし (倍精度) 0.548830

作用反作用の法則を使わないということは、計算量が倍になることを意味します。倍精度の場合、計算量が倍になったとしても作用反作用の法則を使わないやり方(atomicAddを使わない)のほうが早くなります。一方で単精度の場合にはatomicAddがハードウェアサポートされていて、atomicAddを使うオーバーヘッドがそれなりに小さく、作用反作用の法則を使わないで計算量を削減するほうが有利になります。

しかし、あまりにナイーブな実装での比較ですので、ここから二つのカーネルをチューニングしていきながら再度比較していきたいと思います。

チューニングその1 i粒子の運動量書き戻し回数を削減

まず気になるのはi粒子の運動量の書き戻しをループの中で繰り返し行っていることです。これはレジスタに保持しておいてj粒子のループが終わって一回だけグローバルメモリに書き戻すというように変更できます。

作用反作用なし
if (tid < particle_number) {
  const auto qi = q[tid];
  const auto np = number_of_partners[tid];
  const auto kp = pointer[tid];

  Dtype pfx = 0.0, pfy = 0.0, pfz = 0.0; // i粒子の運動量の総和をレジスタに置いておく。
  for (int32_t k = 0; k < np; k++) {
    const auto j = sorted_list[kp + k];
    const auto dx = q[j].x - qi.x;
    const auto dy = q[j].y - qi.y;
    const auto dz = q[j].z - qi.z;
    const auto r2 = dx * dx + dy * dy + dz * dz;
    const auto r6 = r2 * r2 * r2;
    auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
    if (r2 > CL2) df = 0.0;
    pfx += df * dx;
    pfy += df * dy;
    pfz += df * dz;
  }
  p[tid].x += pfx; // 最後に一度だけ書き込みを行う。
  p[tid].y += pfy;
  p[tid].z += pfz;
}
作用反作用あり
const auto tid = threadIdx.x + blockIdx.x * blockDim.x;
if (tid < particle_number) {
  const auto qi = q[tid];
  const auto np = number_of_partners[tid];
  const auto kp = pointer[tid];

  Dtype pfx = 0.0, pfy = 0.0, pfz = 0.0; // i粒子の運動量の総和をレジスタに保持しておく。
    for (int32_t k = 0; k < np; k++) {
      const auto j = sorted_list[kp + k];
      const auto dx = q[j].x - qi.x;
      const auto dy = q[j].y - qi.y;
      const auto dz = q[j].z - qi.z;
      const auto r2 = dx * dx + dy * dy + dz * dz;
      if (r2 > CL2) continue;
      const auto r6 = r2 * r2 * r2;
      const auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
      pfx += df * dx;
      pfy += df * dy;
      pfz += df * dz;
      atomicAdd(&p[j].x, -df * dx);
      atomicAdd(&p[j].y, -df * dy);
      atomicAdd(&p[j].z, -df * dz);
    }
    atomicAdd(&p[tid].x, pfx); // 最後に一度だけ書き戻しを行う。
    atomicAdd(&p[tid].y, pfy);
    atomicAdd(&p[tid].z, pfz);
  }

これでも倍精度の場合をみてみると作用反作用を使わないやり方のほうが速いです。

GPU: Tesla K40t
密度1の時

実行時間 [s]
作用反作用あり 1.345667
作用反作用なし 1.189371

密度0.5の時

実行時間 [s]
作用反作用あり 0.351067
作用反作用なし 0.315363

チューニングその2 近接粒子リストのデータレイアウトの変更

さらにチューニングを進めてみます。
近接粒子リストをメモリからとってくるとき、warp内32スレッドでバラバラのアドレスからとってくるようになっています。コアレスアクセスになっていません。該当箇所のコードを抜粋すると、

...
  const auto qi = q[tid];
  const auto np = number_of_partners[tid];
  const auto kp = pointer[tid];
  Dtype pfx = 0.0, pfy = 0.0, pfz = 0.0;
  for (int32_t k = 0; k < np; k++) { //j粒子についてループ
    const auto j = sorted_list[kp + k]; // warp内32スレッドで連続したメモリにアクセスするわけではない!
...

です。コメントで示したsorted_listという配列へのアクセスはwarp内32スレッドで連続にはなりません。そこで、sorted_listのレイアウトを変更してみます。

どのように変更するのかをわかりやすく図で示すと、変更前は

image

となっています。これを

image

のように変更してコアレスアクセスになるようにします。メモリの消費量が多少大きくなりますが、それさえクリアすれば問題ありません。

結果は

GPU: Tesla K40t
密度1の場合

実行時間 [s]
作用反作用あり 1.331697
作用反作用なし 0.236682

密度0.5の場合

実行時間 [s]
作用反作用あり 0.361904
作用反作用なし 0.087883

これをやるとずいぶん差がついて、作用反作用を使わないやり方がベストなのかなと思えてきます。

チューニングその3 j粒子ループをアンロール

チューニングその2で近接粒子リストのデータレイアウトを変更しましたが、もとのCPUコードをGPGPU移植するさいにはできるだけ避けたいことです。そこで、j粒子ループをアンロールして近接粒子リストのデータレイアウトを変更せずに、コアレスアクセスできるようにしてみました。

何をやるのかというと、j粒子のループを32倍アンロールしてやり、さらに起動するスレッドの数を32倍に増やしてやります。

作用反作用なし
const auto i_ptcl_id = (threadIdx.x + blockIdx.x * blockDim.x) / warpSize;
if (i_ptcl_id < particle_number) {
  const auto lid = lane_id(); // warp内のレーンの番号
  const auto qi = q[i_ptcl_id];
  const auto np = number_of_partners[i_ptcl_id];
  const auto kp = pointer[i_ptcl_id] + lid;
  const int32_t ini_loop = (np / warpSize) * warpSize;

  Vec pf = {0.0};
  if (lid == 0) pf = p[i_ptcl_id];
  int32_t k = 0;
  for (; k < ini_loop; k += warpSize) { // warpSizeでj粒子のループをアンロールする。
    const auto j  = sorted_list[kp + k]; // ここのアクセスはwarp内32スレッドで連続でコアレスアクセス
    const auto dx = q[j].x - qi.x;
    const auto dy = q[j].y - qi.y;
    const auto dz = q[j].z - qi.z;
    const auto r2 = dx * dx + dy * dy + dz * dz;
    const auto r6 = r2 * r2 * r2;
    auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
    if (r2 > CL2) df = 0.0;
    pf.x += df * dx;
    pf.y += df * dy;
    pf.z += df * dz;
  }
...
}
作用反作用あり
const auto i_ptcl_id = (threadIdx.x + blockIdx.x * blockDim.x) / warpSize;
if (i_ptcl_id < particle_number) {
  const auto lid = lane_id();
  const auto qi = q[i_ptcl_id];
  const auto np = number_of_partners[i_ptcl_id];
  const auto kp = pointer[i_ptcl_id] + lid;
  const int32_t ini_loop = (np / warpSize) * warpSize;

  Vec pf = {0.0};
  int32_t k = 0;
  for (; k < ini_loop; k += warpSize) { // warpSizeでj粒子のループをアンロールする。
    const auto j = sorted_list[kp + k]; // ここのアクセスはwarp内32スレッドで連続でコアレスアクセス
    const auto dx = q[j].x - qi.x;
    const auto dy = q[j].y - qi.y;
    const auto dz = q[j].z - qi.z;
    const auto r2 = dx * dx + dy * dy + dz * dz;
    if (r2 <= CL2) {
      const auto r6 = r2 * r2 * r2;
      const auto df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
      pf.x += df * dx;
      pf.y += df * dy;
      pf.z += df * dz;
      atomicAdd(&p[j].x, -df * dx);
      atomicAdd(&p[j].y, -df * dy);
      atomicAdd(&p[j].z, -df * dz);
    }
  }
...
}

j粒子のループをアンロールしてやると、sorted_listへのアクセスはwarp内32スレッドで必ず連続になり、リプレイは発生しなくなります。

最後にi粒子にかかる力の総和を取ってやる必要があります。それはKepler世代で導入されたshuffle命令を用いて実装します。

// warp内32スレッドがレジスタに保持している値を交換し、総和を計算する関数。
template <typename T>                                                           
__device__ __forceinline__ T warp_segment_reduce(T var) {                       
  for (int offset = (warpSize >> 1); offset > 0; offset >>= 1) {                
    var += __shfl_down(var, offset);                                            
  }                                                                             
  return var;                                                                   
}

...

// warp内のpf.xの値の総和を計算する。
pf.x = warp_segment_reduce(pf.x);
pf.y = warp_segment_reduce(pf.y);
pf.z = warp_segment_reduce(pf.z);

// 作用反作用を使う場合にはatomicAddに置き換える必要あり。
if (lid == 0) p[i_ptcl_id] += pf;

この場合の実行時間を比較してやると次のようになります。

Tesla K40tで実行
密度1の場合

実行時間 [s]
作用反作用あり 0.084759
作用反作用なし 0.281980

密度0.5の場合

実行時間 [s]
作用反作用あり 0.039636
作用反作用なし 0.115553

作用反作用を使ったほうが速く、なおかつチューニングその2と比較して倍くらい速いことがわかります。atomicAddのハードウェアサポートがされていない倍精度の場合であっても、作用反作用を使うやり方が速いです。
atomicAddのハードウェアサポートがされているPascal世代のGPUで実行した結果も

Tesla P100で実行
密度1の場合

実行時間 [s]
作用反作用あり 0.021643
作用反作用なし 0.080207

密度0.5の場合

実行時間 [s]
作用反作用あり 0.019731
作用反作用なし 0.041136

となります。「atomic使わない作用反作用使わないやり方がいい」とはもう言えないみたいですね。

ひとまずこれまでのコードとCPUのコードとで速度を比較してみます。比較するコードとしてAVX命令でSIMD化+ソフトウェアパイプライニングしたkaityo256先生のlj_simdを用いました。
CPUはシングルコアでの実行になります。粒子密度は1で粒子数は10万です。

(avxと書いているのがlj_simdの最速コードで実行したときで、それ以外はGPUで実行したときの結果になります。)

Tesla K40tとXeon E5-2680 v3
image

Tesla P100とXeon E5-2680 v3
image

tuning その3の棒グラフがほとんど見えないくらいになりました。Xeon E5-2680 v3 1コアに対して

  • Tesla K40tで実行した場合には35倍の高速化
  • Tesla P100で実行した場合には138倍の高速化

という結果になりました。スレッド並列してNUMA最適化もちゃんと行ったとして、E5-2670 v3は24コアなので、K40tで1.5倍、P100で6倍くらいの差になると思います。(ひとまずlj_simdのマルチコア対応やってみてもう少し真面目な比較をしてみようと思います。)

P100とK40tで4倍違うのはatomicAddがハードウェアサポートされているかどうかと、バンド幅の違いと考えられます。

なぜj粒子ループアンロール+作用反作用ありが最も速いのか、私気になります。

Pascalはともかく、Keplerアーキテクチャでも作用反作用ありのほうが速いというのは結構驚きでした(最初測定ミスかと思いました。)。これについては現在調査中です。

まとめ

LJ力計算に限って言えば、作用反作用無視する利点は無くなってしまったみたいです。Tesla P100の結果をみる限り、ものすごくナイーブに実装したとしても普通にatomicAdd使うほうが高速に計算できるみたいです(単精度、倍精度問わず)。ちゃんと実装して比較してみることが大切です。

追記 (2016/1/26)

OpenMPによるスレッド並列+AVX命令を用いたSIMD化したものを作成して、GPUとの性能比較を行いました。ソースコードはここにあります。
作用反作用を使わないやり方を採用しました(#pragma omp atomicがスカラー変数についてしか使えないのと、#pragma omp criticalが非常に遅くて使い物にならなかったためです)。

実行環境 Intel(R) Xeon(R) CPU E5-2640 v4 @ 2.40GHz 10コア x 2ソケット
実行時間 0.296805 [s]

(スレッド並列せずに1コアで実行した場合の実行時間は4.209338 [sec])

最後にGPU CPUの実行時間を全て比較してみると、

密度1の場合

実行時間 [s] CPUの実行時間からの高速化度合い
Tesla P100 0.021643 14倍
Tesla K40t 0.084759 3.5倍
E5-2640 v4 0.296805 1倍

密度0.5の場合

実行時間 [s] CPUの実行時間からの高速化度合い
Tesla P100 0.019731 4.8倍
Tesla K40t 0.039636 2.4倍
E5-2640 v4 0.094636 1倍

GPUでやると100倍高速化されるというのはさすがに幻想ですが、CPUとGPUをめいいっぱい高速化した上で比較すると、5 ~ 10倍くらいGPUのほうがはやそうです。

各CPU、GPUの理論演算性能、メモリバンド幅、B/F値をリストアップしてみると、

理論演算性能

GFLOPS GB/s B/F値 F/B (B/F値の逆数)
Tesla P100 4036 720 0.18 5.7
Tesla K40t 1430 288 0.20 5.0
E5-2640 v4 768 136.6 0.18 5.6

今回のLJ力計算の演算強度を考えてみると

  • 減算3回 (変位ベクトル計算) $dr_{ij} = r_{i} - r_{j}$
  • 乗算3回 + 足算2回 (距離計算) $dr^{2}_{ij}$
  • 乗算2回 ($dr^{6}$の計算)
  • 乗算2回 ($dr^{14}$の計算)
  • 乗算1回 + 減算1回 (力の分母計算) $24 * dr^{6}_{ij} - 48$
  • 除算1回 $dr^{14}$との割り算
  • 乗算1回 力積のノルム計算 $df * dt$
  • 乗算3回 + 加算3回 (作用反作用を使う場合にはこの2倍の計算) 力積計算 $df * dr_{ij}$

除算を10FLOPとして考えた時、作用反作用を考慮したとき、37FLOPで作用反作用なしで34FLOPになる。

一方で、これらの計算をするのに必要なデータは(i粒子の座標データと運動量データはレジスタにのっているとして、)
* 座標データのロード 3 (4) * 8 Byte = 24 (32) Byte
* 運動量データのロード 2 * 3 (4) * 8 Byte = 48 (64) Byte
となる。(かっこはAVXの_mm256_loadu_pdとかを使って座標を一発でロードする場合に、1成分パディングする場合があるため、それを考慮している。)なお、運動量データのロードは作用反作用を無視する場合には入らない。

よって演算強度は

  • 作用反作用を無視しない場合 37 FLOP / 72 (96) Byte ~ 0.51 (0.38)
  • 作用反作用を無視する場合 34 FLOP / 24 (32) Byte ~ 1.4 (1.1)

になる。ルーフラインモデルで考えると、(ロードするj粒子のデータが全くキャッシュに乗っていないと仮定するのであれば)、メモリバンド幅律速のアプリケーションとなり、バンド幅で性能が決まる。そうすると、二つのGPUで実行した時の実行時間の違いは大まかにはバンド幅の違いとなって表れているはずで、確かにそうなっている。

追記 (2016/1/27)

追記1/26の考察が間違っていると思えてきた。Tesla P100の性能のボトルネックは別にある気がしてきた。
最速カーネルの実行時間をよく調べてみると、doubleとfloatでほとんど変わらない。

Tesla P100で実行
密度1の時

実行時間 [s]
float 0.021128
double 0.021668

ボトルネックがGPUでの計算以外にある気がしてきて、もう少し調べてみると、
ボトルネックがCPU<->GPUのデータ転送にあることがわかった。
今回通信込みで実行時間を計算していたので、通信の時間を見積もってみた。

データ転送の時間をnvvpで調べてみると、通信の時間の合計が0.018sぐらいはあって、実際の計算が0.003sぐらいだと分かった。近接粒子リストのデータが96MBもあるため、転送に相当時間がかかるようだ。
そうすると、実際にGPUで計算している時間で比較すれば、E5-2640 v4のマルチコア実行に比べて100倍ぐらい高速化されている。かなり理想的にキャッシュヒットしているからここまで違いがでると思われる。

metricをnvprofでとってみると、cache hit rateは80%ぐらいだった。

でもCPUも、粒子データが高々5MBくらいだから、ほとんどのデータがキャッシュに乗っていると考えていいような。なぜここまで差が出るのか?