LoginSignup
4
4

More than 5 years have passed since last update.

@kaityo256(ロボ太)

LJの力計算のSIMD化ステップ・バイ・ステップ その5

Last updated at Posted at 2016-11-06

はじめに

せっかくがんばってSIMD化したのに、速度が出なかった。このまま引き下がるわけにはいかない。

コードは
https://github.com/kaityo256/lj_simdstep
においてある。

何が悪かった？

やっぱり、ペア毎にループを回すのはダメだ。今早いコードは「i粒子でソートしたものをSIMD化した」コードだが、スカラで最速のコードは「i粒子でソートしたものをソフトウェアパイプライニングしたもの」なのだから「i粒子でソートしたものをソフトウェアパイプライニングしたものをSIMD化」するしかない。

ただし、ループ長が短いので、前回のようなごまかしは効かない。ちゃんとループの中で先読み計算したものを、ループの外で利用した上で、残りのゴミを計算しないといけない。面倒だがやるしかない。

コード

というわけでやってみた。SIMD化の元となるコードは、「i粒子でソートしたものをソフトウェアパイプライニングしたもの」で、以下のようなコードである。

void
force_sorted_swp(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];
    double pfx = 0;
    double pfy = 0;
    double pfz = 0;
    const int kp = pointer[i];
    int ja = sorted_list[kp];
    double dxa = q[ja][X] - qx_key;
    double dya = q[ja][Y] - qy_key;
    double dza = q[ja][Z] - qz_key;
    double df = 0.0;
    double dxb = 0.0, dyb = 0.0, dzb = 0.0;
    int jb = 0;

    const int np = number_of_partners[i];
    for (int k = kp; k < np + kp; k++) {
      const double dx = dxa;
      const double dy = dya;
      const double dz = dza;
      double r2 = (dx * dx + dy * dy + dz * dz);
      const int j = ja;
      ja = sorted_list[k + 1];
      dxa = q[ja][X] - qx_key;
      dya = q[ja][Y] - qy_key;
      dza = q[ja][Z] - qz_key;
      if (r2 > CL2)continue;
      pfx += df * dxb;
      pfy += df * dyb;
      pfz += df * dzb;
      p[jb][X] -= df * dxb;
      p[jb][Y] -= df * dyb;
      p[jb][Z] -= df * dzb;
      const double r6 = r2 * r2 * r2;
      df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
      jb = j;
      dxb = dx;
      dyb = dy;
      dzb = dz;
    }
    p[jb][X] -= df * dxb;
    p[jb][Y] -= df * dyb;
    p[jb][Z] -= df * dzb;
    p[i][X] += pfx + df * dxb;
    p[i][Y] += pfy + df * dyb;
    p[i][Z] += pfz + df * dzb;
  }
}

これもいい加減長いが、さらに4倍展開してSIMD化する。その作業はこれまでたどった道なので書かない。結果がこちら。

void
force_sorted_swp_intrin(void) {
  const int pn = particle_number;
  const v4df vzero = _mm256_set_pd(0, 0, 0, 0);
  const v4df vcl2 = _mm256_set_pd(CL2, CL2, CL2, CL2);
  const v4df vc24 = _mm256_set_pd(24 * dt, 24 * dt, 24 * dt, 24 * dt);
  const v4df vc48 = _mm256_set_pd(48 * dt, 48 * dt, 48 * dt, 48 * dt);
  for (int i = 0; i < pn; i++) {
    const v4df vqi = _mm256_load_pd((double*)(q + i));
    v4df vpf = _mm256_set_pd(0.0, 0.0, 0.0, 0.0);
    const int kp = pointer[i];
    int ja_1 = sorted_list[kp];
    int ja_2 = sorted_list[kp + 1];
    int ja_3 = sorted_list[kp + 2];
    int ja_4 = sorted_list[kp + 3];
    v4df vqj_1 = _mm256_load_pd((double*)(q + ja_1));
    v4df vdqa_1 = vqj_1 - vqi;
    v4df vqj_2 = _mm256_load_pd((double*)(q + ja_2));
    v4df vdqa_2 = vqj_2 - vqi;
    v4df vqj_3 = _mm256_load_pd((double*)(q + ja_3));
    v4df vdqa_3 = vqj_3 - vqi;
    v4df vqj_4 = _mm256_load_pd((double*)(q + ja_4));
    v4df vdqa_4 = vqj_4 - vqi;

    v4df vdf = _mm256_set_pd(0.0, 0.0, 0.0, 0.0);

    v4df vdqb_1 = _mm256_set_pd(0.0, 0.0, 0.0, 0.0);
    v4df vdqb_2 = _mm256_set_pd(0.0, 0.0, 0.0, 0.0);
    v4df vdqb_3 = _mm256_set_pd(0.0, 0.0, 0.0, 0.0);
    v4df vdqb_4 = _mm256_set_pd(0.0, 0.0, 0.0, 0.0);

    int jb_1 = 0, jb_2 = 0, jb_3 = 0, jb_4 = 0;
    const int np = number_of_partners[i];
    for (int k = 0; k < (np / 4) * 4; k += 4) {
      const int j_1 = ja_1;
      const int j_2 = ja_2;
      const int j_3 = ja_3;
      const int j_4 = ja_4;
      v4df vdq_1 = vdqa_1;
      v4df vdq_2 = vdqa_2;
      v4df vdq_3 = vdqa_3;
      v4df vdq_4 = vdqa_4;

      ja_1 = sorted_list[kp + k + 4];
      ja_2 = sorted_list[kp + k + 5];
      ja_3 = sorted_list[kp + k + 6];
      ja_4 = sorted_list[kp + k + 7];

      v4df vr2s_1 = vdq_1 * vdq_1;
      v4df vr2t_1 = _mm256_permute4x64_pd(vr2s_1, 201);
      v4df vr2u_1 = _mm256_permute4x64_pd(vr2s_1, 210);
      v4df vr2_1 = vr2s_1 + vr2t_1 + vr2u_1;

      v4df vr2s_2 = vdq_2 * vdq_2;
      v4df vr2t_2 = _mm256_permute4x64_pd(vr2s_2, 201);
      v4df vr2u_2 = _mm256_permute4x64_pd(vr2s_2, 210);
      v4df vr2_2 = vr2s_2 + vr2t_2 + vr2u_2;

      v4df vr2s_3 = vdq_3 * vdq_3;
      v4df vr2t_3 = _mm256_permute4x64_pd(vr2s_3, 201);
      v4df vr2u_3 = _mm256_permute4x64_pd(vr2s_3, 210);
      v4df vr2_3 = vr2s_3 + vr2t_3 + vr2u_3;

      v4df vr2s_4 = vdq_4 * vdq_4;
      v4df vr2t_4 = _mm256_permute4x64_pd(vr2s_4, 201);
      v4df vr2u_4 = _mm256_permute4x64_pd(vr2s_4, 210);
      v4df vr2_4 = vr2s_4 + vr2t_4 + vr2u_4;

      v4df vdf_1 = _mm256_permute4x64_pd(vdf, 0);
      v4df vdf_2 = _mm256_permute4x64_pd(vdf, 85);
      v4df vdf_3 = _mm256_permute4x64_pd(vdf, 170);
      v4df vdf_4 = _mm256_permute4x64_pd(vdf, 255);

      vqj_1 = _mm256_load_pd((double*)(q + ja_1));
      vdqa_1 = vqj_1 - vqi;
      vpf += vdf_1 * vdqb_1;

      v4df vpjb_1 = _mm256_load_pd((double*)(p + jb_1));
      vpjb_1 -= vdf_1 * vdqb_1;
      _mm256_store_pd((double*)(p + jb_1), vpjb_1);

      vqj_2 = _mm256_load_pd((double*)(q + ja_2));
      vdqa_2 = vqj_2 - vqi;
      vpf += vdf_2 * vdqb_2;

      v4df vpjb_2 = _mm256_load_pd((double*)(p + jb_2));
      vpjb_2 -= vdf_2 * vdqb_2;
      _mm256_store_pd((double*)(p + jb_2), vpjb_2);

      vqj_3 = _mm256_load_pd((double*)(q + ja_3));
      vdqa_3 = vqj_3 - vqi;
      vpf += vdf_3 * vdqb_3;

      v4df vpjb_3 = _mm256_load_pd((double*)(p + jb_3));
      vpjb_3 -= vdf_3 * vdqb_3;
      _mm256_store_pd((double*)(p + jb_3), vpjb_3);

      vqj_4 = _mm256_load_pd((double*)(q + ja_4));
      vdqa_4 = vqj_4 - vqi;
      vpf += vdf_4 * vdqb_4;

      v4df vpjb_4 = _mm256_load_pd((double*)(p + jb_4));
      vpjb_4 -= vdf_4 * vdqb_4;
      _mm256_store_pd((double*)(p + jb_4), vpjb_4);

      v4df vr2_13 = _mm256_unpacklo_pd(vr2_1, vr2_3);
      v4df vr2_24 = _mm256_unpacklo_pd(vr2_2, vr2_4);
      v4df vr2 = _mm256_shuffle_pd(vr2_13, vr2_24, 12);
      v4df vr6 = vr2 * vr2 * vr2;
      vdf = (vc24 * vr6 - vc48) / (vr6 * vr6 * vr2);
      v4df mask = vcl2 - vr2;
      vdf = _mm256_blendv_pd(vdf, vzero, mask);

      jb_1 = j_1;
      jb_2 = j_2;
      jb_3 = j_3;
      jb_4 = j_4;
      vdqb_1 = vdq_1;
      vdqb_2 = vdq_2;
      vdqb_3 = vdq_3;
      vdqb_4 = vdq_4;
    }
    v4df vdf_1 = _mm256_permute4x64_pd(vdf, 0);
    v4df vdf_2 = _mm256_permute4x64_pd(vdf, 85);
    v4df vdf_3 = _mm256_permute4x64_pd(vdf, 170);
    v4df vdf_4 = _mm256_permute4x64_pd(vdf, 255);

    v4df vpjb_1 = _mm256_load_pd((double*)(p + jb_1));
    vpjb_1 -= vdf_1 * vdqb_1;
    _mm256_store_pd((double*)(p + jb_1), vpjb_1);

    v4df vpjb_2 = _mm256_load_pd((double*)(p + jb_2));
    vpjb_2 -= vdf_2 * vdqb_2;
    _mm256_store_pd((double*)(p + jb_2), vpjb_2);

    v4df vpjb_3 = _mm256_load_pd((double*)(p + jb_3));
    vpjb_3 -= vdf_3 * vdqb_3;
    _mm256_store_pd((double*)(p + jb_3), vpjb_3);

    v4df vpjb_4 = _mm256_load_pd((double*)(p + jb_4));
    vpjb_4 -= vdf_4 * vdqb_4;
    _mm256_store_pd((double*)(p + jb_4), vpjb_4);

    v4df vpi = _mm256_load_pd((double*)(p + i));
    vpf += vdf_1 * vdqb_1;
    vpf += vdf_2 * vdqb_2;
    vpf += vdf_3 * vdqb_3;
    vpf += vdf_4 * vdqb_4;
    vpi += vpf;
    _mm256_store_pd((double*)(p + i), vpi);
    const double qix = q[i][X];
    const double qiy = q[i][Y];
    const double qiz = q[i][Z];
    double pfx = 0.0;
    double pfy = 0.0;
    double pfz = 0.0;
    for (int k = (np / 4) * 4; k < np; k++) {
      const int j = sorted_list[k + kp];
      double dx = q[j][X] - qix;
      double dy = q[j][Y] - qiy;
      double dz = q[j][Z] - qiz;
      double r2 = (dx * dx + dy * dy + dz * dz);
      double r6 = r2 * r2 * r2;
      double 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[j][X] -= df * dx;
      p[j][Y] -= df * dy;
      p[j][Z] -= df * dz;
    }
    p[i][X] += pfx;
    p[i][Y] += pfy;
    p[i][Z] += pfz;
  }
}

うん、ちょっと気が滅入るコード量。普通の気分なら絶対やらなかったな。

結果

というわけで結果。計算しなおしたから、これまでと微妙に数値が変わるけど気にしないこと。

計算手法 実行時間[sec] 備考
pair 8.175502 ペアごとにループをまわしたもの
pair_swp 6.994964 上記をSWPしたもの
pair_swp_intrin 4.273858 上記を4倍展開してSIMD化したもの
sorted 7.072854 i粒子でソートしたもの
sorted_intrin 3.880328 i粒子でソートしたものをSIMD化したもの
sorted_swp 7.072854 i粒子でソートしたものをSWPしたもの
sorted_swp_intrin 3.278851 i粒子でソートしたものをSWPしたものをSIMD化したもの

まぁ若干早くなりましたな。上記は高密度(数密度1.0)の結果だが、低密度(数密度0.5)の結果はこんな感じ。

計算手法 実行時間[sec] 備考
pair 1.692347 ペアごとにループをまわしたもの
pair_swp 1.747897 上記をSWPしたもの
pair_swp_intrin 1.197281 上記を4倍展開してSIMD化したもの
sorted 1.523837 i粒子でソートしたもの
sorted_intrin 1.232721 i粒子でソートしたものをSIMD化したもの
sorted_swp 1.069575 i粒子でソートしたものをSWPしたもの
sorted_swp_intrin 0.965486 i粒子でソートしたものをSWPしたものをSIMD化したもの

まぁ、少しだけSWP(ソフトウェアパイプライニング)の効果はあったかな。

コードは以下に。

まとめ

というわけでSIMD化してみた。少なくとも馬鹿SIMD化は、データ構造とループ構造が決まっていれば単純作業なので難しくはない。ただし、SIMD化後の性能向上のかなりの部分は「SIMD化前のループ構造」で決まることがわかった。あ、

SIMD化の性能はSIMD化前に決まっている

ってなんか名言っぽくないですか？

おしまい・・・のはずが続いた

4
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
4