AVX命令を用いたVerlet list構築のSIMD化 その4

はじめに

前回はループ展開する場所の変更により性能を向上させた。これ以上の高速化を行うことはできないか?

Verlet list構築のSIMD化シリーズ
* その1 最内ループ4倍展開+馬鹿SIMD+逐次ストア
* その2 逐次ストアから64~256bitストアに変更
* その3 展開するループを変更し、メモリアクセス最適化
* その4 最終回

SIMD化の割合を高める

現在のループ構造は

for (すべてのセルについて) {
  for (int i = 0; i < (iloop / 4) * 4; i += 4) {
    v4df vqia = ...;
    v4df vqib = ...;
    v4df vqic = ...;
    v4df vqid = ...;
    for (int j = 0; j < jloop; j++) {
      // Verlet list構築
      v4df vqj = ...;
      //64bit ~ 256bitの粒子ペアデータをストア
    }
  }

  // i粒子ループを展開することによって生じるゴミループ
  for (int i = (iloop / 4) * 4; i < iloop; i++) {
    for (int j = 0; j < jloop; j++) { // この部分はSIMD化できる!
      // Verlet list構築
    }
  }
}

iループの展開に際して発生するゴミループは結構重い。jloop x (iloop % 3)分のループが余っている。
このあまりの部分についてはj loopについてSIMD化できる。

4x4行列の転置の回数を減らす

これまで

const auto ri_a = q[i_a];
const auto ri_b = q[i_b];
const auto ri_c = q[i_c];
const auto ri_d = q[i_d];
for (int j = 0; j < jloop; j++) { // 最内ループ
  const auto drij_a = ri_a - q[j];
  const auto drij_b = ri_b - q[j];
  const auto drij_c = ri_c - q[j];
  const auto drij_d = ri_d - q[j];

  ...

  // 4x4 行列の転置を行う。
  /* drij_a = {dxa, dya, dza, 0}  4x4transpose dx_abcd = {dxa, dxb, dxc, dxd}
     drij_b = {dxb, dyb, dzb, 0}        ->     dy_abcd = {dya, dyb, dyc, dyd}
     drij_c = {dxc, dyc, dzc, 0}               dz_abcd = {dza, dzb, dzc, dzd}
     drij_d = {dxd, dyd, dzd, 0}               dw_abcd = {  0,   0,   0,   0} */
  //

  ...

  // 変位ベクトルの二乗
  const auto drij2_abcd = dx_abcd * dx_abcd + dy_abcd * dy_abcd + dz_abcd * dz_abcd;

  // Verlet listをストア
  // drij2_abcdを計算したあとで、変位ベクトルdrij_{abcd}のデータは必要ない。

  ...
}

のように変位ベクトルを計算したのち、行列の転置を行いそこから距離の二乗を計算するようにしていた。
しかし、力の計算の場合と違い、変位ベクトルdrij_a drij_b drij_c drij_dは、距離の二乗の計算が終わってしまえば、最内ループの以降の処理に必要ない。上記のコメントに示しているdx_abcd dy_abcd dz_abcdを直接計算することができれば距離の二乗の計算はできる。

そこで、以下の疑似コードに示すように

const auto ri_a = q[i_a];
const auto ri_b = q[i_b];
const auto ri_c = q[i_c];
const auto ri_d = q[i_d];

...

// 最初にri_a ri_b ri_c ri_dの 4x4 行列の転置を行う。
/* ri_a = {xia, yia, zia, 0}  4x4transpose xi_abcd = {xia, xib, xic, xid}
   ri_b = {xib, yib, zib, 0}        ->     yi_abcd = {yia, yib, yic, yid}
   ri_c = {xic, yic, zic, 0}               zi_abcd = {zia, zib, zic, zid}
   ri_d = {xid, yid, zid, 0}               wi_abcd = {  0,   0,   0,   0} */
//

...

for (int j = 0; j < jloop; j++) { // 最内ループ
  v4df xj = _mm256_set1_pd(q[j].x); // xj = {q[j].x, q[j].x, q[j].x, q[j].x}
  v4df yj = _mm256_set1_pd(q[j].y); // yj = {q[j].y, q[j].y, q[j].y, q[j].y}
  v4df zj = _mm256_set1_pd(q[j].z); // zj = {q[j].z, q[j].z, q[j].z, q[j].z}

  // drij_a drij_b drij_c drij_dを作らない。
  // 変位ベクトルの各成分(x, y, z)が格納されたv4df型のデータを用意する。
  const auto dx_abcd = xi_abcd - xj;
  const auto dy_abcd = yi_abcd - yj;
  const auto dz_abcd = zi_abcd - zj;

  const auto drij2_abcd = dx_abcd * dx_abcd + dy_abcd * dy_abcd + dz_abcd * dz_abcd;

  // Verlet listをストア

  ...
}

に変更する。アンロールして出てきたi粒子の4x256bitの座標データを最初に一回だけ転置してx{yz}i_abcdのデータを作成したところが大きな変更点である。これによりdrij_a drij_b drij_c drij_dの計算を経由せず、直接距離の二乗の計算を行う。

結果

ソースコードはここのMakePairListFusedLoopSIMD4x1AllTransである。

実行環境は
* CPU: Intel(R) Xeon(R) CPU E5-2680 v3 @ 2.50GHz
* コンパイラ: icpc (ICC) 16.0.4 20160811
(前回実行環境書くの忘れていました。すみません。)

100回近接粒子リストを構築するのに要した時間
(ただし、密度1で粒子数10万での結果)

実行時間
reference 22593[ms]
最内ループを4倍展開し距離計算をSIMD化して逐次ストア 23948[ms]
上記に加えストアをSIMD化 17794[ms]
最内ループより一つ外側を4倍展開して距離計算とストアをSIMD化 15629[ms]
上記に加えSIMD化率引き上げ、シャッフル命令を削減したもの 12173[ms]
上記にソフトウェアパイプライニングを加えたもの 12002[ms]

目標としていた2倍弱くらいにはなった。

補足 vgather系命令を用いた場合

const auto ri_a = q[i_a];
const auto ri_b = q[i_b];
const auto ri_c = q[i_c];
const auto ri_d = q[i_d];

...

// 最初にri_a ri_b ri_c ri_dの 4x4 行列の転置を行う。
/* ri_a = {xia, yia, zia, 0}  4x4transpose xi_abcd = {xia, xib, xic, xid}
   ri_b = {xib, yib, zib, 0}        ->     yi_abcd = {yia, yib, yic, yid}
   ri_c = {xic, yic, zic, 0}               zi_abcd = {zia, zib, zic, zid}
   ri_d = {xid, yid, zid, 0}               wi_abcd = {  0,   0,   0,   0} */
//

の部分で、xi_abcd yi_abcd zi_abcdを得るのに転置を使わずにvgather系の命令で得ることもできる。例えば、

v4si vindex = _mm_set_epi32(4 * i_d, 4 * i_c, 4 * i_b, 4 * i_a);
v4df xi_abcd = _mm256_i32gather_pd(&(q[0].x), vindex, 8);
v4df yi_abcd = _mm256_i32gather_pd(&(q[0].y), vindex, 8);
v4df zi_abcd = _mm256_i32gather_pd(&(q[0].z), vindex, 8);

のようにして。

しかし、簡単な例で試した限りではvgather使わない方が20%ほど早かったので、行列を転置するやり方を採用した。