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

# はじめに

LJのSIMD化の試み。今回はソフトウェアパイプライニングしたコードを４倍展開する。

* [その0](http://qiita.com/kaityo256/items/5be3ce71a6ff9522a0e6) 
* [その1](http://qiita.com/kaityo256/items/4d34bf144122b78fcf49)
* [その2](http://qiita.com/kaityo256/items/e6c41b94c5cc5189c5cc) 
* [その3](http://qiita.com/kaityo256/items/7f16ab7f2281bb03e83e) ←イマココ
* [その3.5](http://qiita.com/kaityo256/items/b145b2d83becccd0aa53)
* [その4](http://qiita.com/kaityo256/items/636c5a407f5a6c61bb21)
* [その5](http://qiita.com/kaityo256/items/dacbe5a16abbc83eb107)
* [その6](http://qiita.com/kaityo256/items/34c371c6d040be1bdc30)

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

※ っていうかタグにAVX2とか書いてあるのに、ここまでAVX2命令をちっとも使ってないのが詐欺っぽい。

# 4倍展開

ループをソフトウェアパイプライニングできたので、これを馬鹿SIMD化のために4倍展開する。ただし、ソフトウェアパイプライニングしているので端の処理面倒になる。真面目にやるならループ内で計算した結果をループ外でちゃんと使ってから、ゴミの処理をしなければならないが、面倒なので一番端は先読みした結果を捨てて、再度計算してしまおう。どうせループ長がすごく長い(今回のケースで7839886回転)ので、そのコストは無視できる。

そういった方針で作った4倍展開のコードはこちら。これからSIMD化するため、名前をforce_pair_swp_intrinとしてある。

```cpp
void
force_pair_swp_intrin(void){
  int k = 0;
  int i_a1 = i_particles[k];
  int j_a1 = j_particles[k];
  int i_a2 = i_particles[k+1];
  int j_a2 = j_particles[k+1];
  int i_a3 = i_particles[k+2];
  int j_a3 = j_particles[k+2];
  int i_a4 = i_particles[k+3];
  int j_a4 = j_particles[k+3];
  double dx_b1 = q[j_a1][X] - q[i_a1][X];
  double dy_b1 = q[j_a1][Y] - q[i_a1][Y];
  double dz_b1 = q[j_a1][Z] - q[i_a1][Z];
  double dx_b2 = q[j_a2][X] - q[i_a2][X];
  double dy_b2 = q[j_a2][Y] - q[i_a2][Y];
  double dz_b2 = q[j_a2][Z] - q[i_a2][Z];
  double dx_b3 = q[j_a3][X] - q[i_a3][X];
  double dy_b3 = q[j_a3][Y] - q[i_a3][Y];
  double dz_b3 = q[j_a3][Z] - q[i_a3][Z];
  double dx_b4 = q[j_a4][X] - q[i_a4][X];
  double dy_b4 = q[j_a4][Y] - q[i_a4][Y];
  double dz_b4 = q[j_a4][Z] - q[i_a4][Z];
  double dx_a1,dy_a1,dz_a1;
  double dx_a2,dy_a2,dz_a2;
  double dx_a3,dy_a3,dz_a3;
  double dx_a4,dy_a4,dz_a4;
  int i_b1, j_b1;
  int i_b2, j_b2;
  int i_b3, j_b3;
  int i_b4, j_b4;
  double df_1,df_2,df_3,df_4;
  for(k=4;k<(number_of_pairs)/4*4;k+=4){
    dx_a1 = dx_b1; 
    dy_a1 = dy_b1; 
    dz_a1 = dz_b1; 
    dx_a2 = dx_b2; 
    dy_a2 = dy_b2; 
    dz_a2 = dz_b2; 
    dx_a3 = dx_b3; 
    dy_a3 = dy_b3; 
    dz_a3 = dz_b3; 
    dx_a4 = dx_b4; 
    dy_a4 = dy_b4; 
    dz_a4 = dz_b4; 

    i_b1 = i_particles[k];
    j_b1 = j_particles[k];
    i_b2 = i_particles[k+1];
    j_b2 = j_particles[k+1];
    i_b3 = i_particles[k+2];
    j_b3 = j_particles[k+2];
    i_b4 = i_particles[k+3];
    j_b4 = j_particles[k+3];

    dx_b1 = q[j_b1][X] - q[i_b1][X];
    dy_b1 = q[j_b1][Y] - q[i_b1][Y];
    dz_b1 = q[j_b1][Z] - q[i_b1][Z];
    dx_b2 = q[j_b2][X] - q[i_b2][X];
    dy_b2 = q[j_b2][Y] - q[i_b2][Y];
    dz_b2 = q[j_b2][Z] - q[i_b2][Z];
    dx_b3 = q[j_b3][X] - q[i_b3][X];
    dy_b3 = q[j_b3][Y] - q[i_b3][Y];
    dz_b3 = q[j_b3][Z] - q[i_b3][Z];
    dx_b4 = q[j_b4][X] - q[i_b4][X];
    dy_b4 = q[j_b4][Y] - q[i_b4][Y];
    dz_b4 = q[j_b4][Z] - q[i_b4][Z];

    const double r2_1 = (dx_a1 * dx_a1 + dy_a1 * dy_a1 + dz_a1 * dz_a1);
    const double r2_2 = (dx_a2 * dx_a2 + dy_a2 * dy_a2 + dz_a2 * dz_a2);
    const double r2_3 = (dx_a3 * dx_a3 + dy_a3 * dy_a3 + dz_a3 * dz_a3);
    const double r2_4 = (dx_a4 * dx_a4 + dy_a4 * dy_a4 + dz_a4 * dz_a4);
    const double r6_1 = r2_1 * r2_1 * r2_1;
    const double r6_2 = r2_2 * r2_2 * r2_2;
    const double r6_3 = r2_3 * r2_3 * r2_3;
    const double r6_4 = r2_4 * r2_4 * r2_4;
    df_1 = ((24.0 * r6_1 - 48.0) / (r6_1 * r6_1 * r2_1)) * dt;
    df_2 = ((24.0 * r6_2 - 48.0) / (r6_2 * r6_2 * r2_2)) * dt;
    df_3 = ((24.0 * r6_3 - 48.0) / (r6_3 * r6_3 * r2_3)) * dt;
    df_4 = ((24.0 * r6_4 - 48.0) / (r6_4 * r6_4 * r2_4)) * dt;

    if (r2_1 > CL2) df_1=0.0;
    if (r2_2 > CL2) df_2=0.0;
    if (r2_3 > CL2) df_3=0.0;
    if (r2_4 > CL2) df_4=0.0;
    p[i_a1][X] += df_1 * dx_a1;
    p[i_a1][Y] += df_1 * dy_a1;
    p[i_a1][Z] += df_1 * dz_a1;
    p[j_a1][X] -= df_1 * dx_a1;
    p[j_a1][Y] -= df_1 * dy_a1;
    p[j_a1][Z] -= df_1 * dz_a1;

    p[i_a2][X] += df_2 * dx_a2;
    p[i_a2][Y] += df_2 * dy_a2;
    p[i_a2][Z] += df_2 * dz_a2;
    p[j_a2][X] -= df_2 * dx_a2;
    p[j_a2][Y] -= df_2 * dy_a2;
    p[j_a2][Z] -= df_2 * dz_a2;

    p[i_a3][X] += df_3 * dx_a3;
    p[i_a3][Y] += df_3 * dy_a3;
    p[i_a3][Z] += df_3 * dz_a3;
    p[j_a3][X] -= df_3 * dx_a3;
    p[j_a3][Y] -= df_3 * dy_a3;
    p[j_a3][Z] -= df_3 * dz_a3;

    p[i_a4][X] += df_4 * dx_a4;
    p[i_a4][Y] += df_4 * dy_a4;
    p[i_a4][Z] += df_4 * dz_a4;
    p[j_a4][X] -= df_4 * dx_a4;
    p[j_a4][Y] -= df_4 * dy_a4;
    p[j_a4][Z] -= df_4 * dz_a4;

    i_a1 = i_b1;
    j_a1 = j_b1;
    i_a2 = i_b2;
    j_a2 = j_b2;
    i_a3 = i_b3;
    j_a3 = j_b3;
    i_a4 = i_b4;
    j_a4 = j_b4;
  }
  for(k=(number_of_pairs)/4*4-4;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;
  }
}
```

アホみたいに長くなりますね。これ、後で変数の付け方が悪いことがわかったので、後で4倍展開する時には別の命名法を採用します(予告)。端の処理を手抜きしたのが、最後の

```cpp
  for(k=(number_of_pairs)/4*4-4;k<number_of_pairs;k++){//★
```
にあらわれている。本来は`-4`は不要なのだが、先読みした結果を捨てて再計算している。先に述べたように、これは今回は問題にはならない。

# 結果

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

なんか実行時間がやけに長くなってるけど、気にしないでこれをSIMD化することにしよう。

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

[次回へ続く](http://qiita.com/kaityo256/items/b145b2d83becccd0aa53)。