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

はじめに

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

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

コード

ペア毎の力計算コードは、ループが一重ループになっている。コードはこんな感じ。

void
force_pair(void){
  for(int k=0;k<number_of_pairs;k++){
    const int i = i_particles[k];
    const 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);
    if (r2 > CL2) continue;
    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;
  }
}

極めてシンプル。kがループインデックスで、それがペア番号になっている。それからi粒子の番号とj粒子の番号を取ってきて、相対座標を計算し、距離を計算し、力を計算し、力積を書き戻す、という手順。ただし、Bookkeeping法という手法を用いているので、相互作用範囲外にあるペアもペアリストに登録されている。それらは距離を調べて弾いている。

これをソフトウェアパイプライニングしよう。いろいろ試行錯誤したが、「距離の計算の直前」(★のあるところ)で切って、その前後を重ねるのが一番良さそうである。そのように修正したコードは以下の通り。

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;
}

ただ半分ずらしただけで、ループの前後に出てくる「はみ出し」の処理が面倒なこと以外は特に難しいことはない。

結果

make testは、力を計算し、最初の5粒子と最後の5粒子の運動量を表示して、他の手法と結果が一致するか確認する。これで結果が正しいことをまず確認(基本)。

そして、肝心の計算時間だが・・・

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

遅くなってしまった。計算は独立なのだが、continue文があるのと、インデックスを共通のものを使っているので、独立に計算できることをコンパイラが判断できなかったようだ。

今回のコードはこちら。 https://github.com/kaityo256/lj_simdstep/blob/master/step1/force.cpp

次回へ続く。