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

はじめに

LJのSIMD化の試み。まずはペアのループをソフトウェアパイプライニングする。

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

前回までのあらすじ

ソフトウェアパイプライニングしたは良いけど遅くなってしまった。おそらくcontinue文が悪いのであろう。また、インデックスを共有しているのも良くない気がする。そこで、continueをやめてゼロクリアし、さらに独立に計算できる力の計算を明示的にループの前半部分に移動してみる。インデックスもi,jの再利用をやめ、i_a, j_a, i_b, j_bと明示的に別の場所であることを明示する。

修正前はこんな感じ。

void
force_pair_swp(void){
  int k = 0;
  int i = i_particles[k];
  int j = j_particles[k];
  double dx = q[j][X] - q[i][X];
  double dy = q[j][Y] - q[i][Y];
  double dz = q[j][Z] - q[i][Z];
  double r2 = (dx * dx + dy * dy + dz * dz);
  double r6, df;
  for(k=1;k<number_of_pairs;k++){
    r2 = (dx * dx + dy * dy + dz * dz);
    r6 = r2 * r2 * r2;
    df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
    if (r2 > CL2) df=0.0;
    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;
    i = i_particles[k];
    j = j_particles[k];
    dx = q[j][X] - q[i][X];
    dy = q[j][Y] - q[i][Y];
    dz = q[j][Z] - q[i][Z];
  }
  r2 = (dx * dx + dy * dy + dz * dz);
  r6 = r2 * r2 * r2;
  df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
  if (r2 > CL2) df=0.0;
  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;
}

以下が修正後。

void
force_pair_swp(void){
  int k = 0;
  int i_a = i_particles[k];
  int j_a = j_particles[k];
  double dx_b = q[j_a][X] - q[i_a][X];
  double dy_b = q[j_a][Y] - q[i_a][Y];
  double dz_b = q[j_a][Z] - q[i_a][Z];
  double dx_a,dy_a,dz_a;
  int i_b, j_b;
  double df;
  for(k=1;k<number_of_pairs;k++){
    dx_a = dx_b; 
    dy_a = dy_b; 
    dz_a = dz_b; 
    i_b = i_particles[k];
    j_b = j_particles[k];
    dx_b = q[j_b][X] - q[i_b][X];
    dy_b = q[j_b][Y] - q[i_b][Y];
    dz_b = q[j_b][Z] - q[i_b][Z];
    const double r2 = (dx_a * dx_a + dy_a * dy_a + dz_a * dz_a);
    const double r6 = r2 * r2 * r2;
    df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
    if (r2 > CL2) df=0.0;
    p[i_a][X] += df * dx_a;
    p[i_a][Y] += df * dy_a;
    p[i_a][Z] += df * dz_a;
    p[j_a][X] -= df * dx_a;
    p[j_a][Y] -= df * dy_a;
    p[j_a][Z] -= df * dz_a;
    i_a = i_b;
    j_a = j_b;
  }
  dx_a = dx_b; 
  dy_a = dy_b; 
  dz_a = dz_b; 
  const double r2 = (dx_a * dx_a + dy_a * dy_a + dz_a * dz_a);
  const double r6 = r2 * r2 * r2;
  df = ((24.0 * r6 - 48.0) / (r6 * r6 * r2)) * dt;
  if (r2 > CL2) df=0.0;
  p[i_a][X] += df * dx_a;
  p[i_a][Y] += df * dy_a;
  p[i_a][Z] += df * dz_a;
  p[j_a][X] -= df * dx_a;
  p[j_a][Y] -= df * dy_a;
  p[j_a][Z] -= df * dz_a;
}

結果

計算手法 実行時間[sec] 備考
pair 8.176446 ペアごとにループをまわしたもの
pair_swp 6.995374 上記をソフトウェアパイプライニングしたもの

無事に早くなった。perf statで調べてみると、pairのIPCが1.33、pair_swpのIPCが1.99 と向上しており、それが命令数の増加を上回っているため早くなったようだ。これをSIMD化することにしよう。

ここまでのコードは以下においてある。 https://github.com/kaityo256/lj_simdstep/tree/master/step2

次回へ続く。