3
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Last updated at Posted at 2016-11-06

はじめに

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

次回へ続く

3
3
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
3
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?