gccとclangの浮動小数点演算(SSE2)に関する最適化具合の違い

先日『[gccとclangのラムダ式に関する最適化具合の違い](http://qiita.com/aokomoriuta/items/517bf29f0cb62940fb31)』というのを見つけたばかりですが（結果的にラムダ式というよりは参照に関する最適化に近かったですが）、今度は顕著に実行時間が1.5倍以上違う（clangの方が1.5倍速い）結果が出たので置いておきます。

実行可能なコードは[bitbucket](https://bitbucket.org/aokomoriuta/simd/src/simple/cpp/simd/simd.cpp)にあるので試してみてください（普通の``make``だとgcc、``make -f Makefile_clang``でclangになります）。
私の以下の環境だと、gccが約5.1[s]で、clangが約3.4[s]でした。

* Ubuntu 15.04 (GNU/Linux 3.19.0-26-generic x86_64)
* Intel(R) Core(TM) i7-3770K CPU @ 3.50GHz
* DDR3-1600 4GBx2
* gcc (Ubuntu 4.9.2-10ubuntu13) 4.9.2
* clang version 3.6.0-2ubuntu1 (tags/RELEASE_360/final) (based on LLVM 3.6.0)

これだとよく分からなかったので色々と削っていくと、最終的に以下の様なコード

```force.cpp
#include <cmath>

struct Vector4
{
	double data[4];
};

void Force(const Vector4 x[], Vector4 f[], const std::size_t n)
{
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			if (i != j)
			{
				const double x0 = x[i].data[0];
				const double x1 = x[i].data[1];
				const double x2 = x[i].data[2];
				const double x3 = x[i].data[3];
				const double r2 = (x0*x0 + x1*x1) + (x2*x2 + x3*x3);

				const double f0 = f[i].data[0] + x0/r2;
				const double f1 = f[i].data[1] + x1/r2;
				const double f2 = f[i].data[2] + x2/r2;
				const double f3 = f[i].data[3] + x3/r2;
				f[i].data[0] = f0;
				f[i].data[1] = f1;
				f[i].data[2] = f2;
				f[i].data[3] = f3;
			}
		}
	}
}

```

を、

- ``g++ -c -O3 force.cpp && objdump -M intel -S force.o``
- ``clang++ -c -O3 force.cpp && objdump -M intel -S force.o``

でコンパイルして中身を覗いてみると、


```cpp-objdump:gcc.txt
   0:   48 85 d2                test   rdx,rdx
   3:   0f 84 a2 00 00 00       je     ab <_Z6NormalPK7Vector4PS_m+0xab>
   9:   31 c9                   xor    ecx,ecx
   b:   0f 1f 44 00 00          nop    DWORD PTR [rax+rax*1+0x0]
  10:   31 c0                   xor    eax,eax
  12:   66 0f 1f 44 00 00       nop    WORD PTR [rax+rax*1+0x0]
  18:   39 c1                   cmp    ecx,eax
  1a:   74 75                   je     91 <_Z6NormalPK7Vector4PS_m+0x91>
  1c:   f2 0f 10 2f             movsd  xmm5,QWORD PTR [rdi]
  20:   f2 0f 10 67 08          movsd  xmm4,QWORD PTR [rdi+0x8]
  25:   66 0f 28 f5             movapd xmm6,xmm5
  29:   f2 0f 10 5f 10          movsd  xmm3,QWORD PTR [rdi+0x10]
  2e:   66 0f 28 c4             movapd xmm0,xmm4
  32:   f2 0f 59 f5             mulsd  xmm6,xmm5
  36:   f2 0f 59 c4             mulsd  xmm0,xmm4
  3a:   f2 0f 10 57 18          movsd  xmm2,QWORD PTR [rdi+0x18]
  3f:   66 0f 28 ca             movapd xmm1,xmm2
  43:   f2 0f 58 f0             addsd  xmm6,xmm0
  47:   66 0f 28 c3             movapd xmm0,xmm3
  4b:   f2 0f 59 ca             mulsd  xmm1,xmm2
  4f:   f2 0f 59 c3             mulsd  xmm0,xmm3
  53:   f2 0f 58 c1             addsd  xmm0,xmm1
  57:   f2 0f 58 c6             addsd  xmm0,xmm6
  5b:   f2 0f 5e e0             divsd  xmm4,xmm0
  5f:   f2 0f 5e d8             divsd  xmm3,xmm0
  63:   f2 0f 58 66 08          addsd  xmm4,QWORD PTR [rsi+0x8]
  68:   f2 0f 11 66 08          movsd  QWORD PTR [rsi+0x8],xmm4
  6d:   f2 0f 5e d0             divsd  xmm2,xmm0
  71:   f2 0f 58 5e 10          addsd  xmm3,QWORD PTR [rsi+0x10]
  76:   f2 0f 11 5e 10          movsd  QWORD PTR [rsi+0x10],xmm3
  7b:   f2 0f 5e e8             divsd  xmm5,xmm0
  7f:   f2 0f 58 56 18          addsd  xmm2,QWORD PTR [rsi+0x18]
  84:   f2 0f 11 56 18          movsd  QWORD PTR [rsi+0x18],xmm2
  89:   f2 0f 58 2e             addsd  xmm5,QWORD PTR [rsi]
  8d:   f2 0f 11 2e             movsd  QWORD PTR [rsi],xmm5
  91:   83 c0 01                add    eax,0x1
  94:   39 d0                   cmp    eax,edx
  96:   75 80                   jne    18 <_Z6NormalPK7Vector4PS_m+0x18>
  98:   83 c1 01                add    ecx,0x1
  9b:   48 83 c7 20             add    rdi,0x20
  9f:   48 83 c6 20             add    rsi,0x20
  a3:   39 d1                   cmp    ecx,edx
  a5:   0f 85 65 ff ff ff       jne    10 <_Z6NormalPK7Vector4PS_m+0x10>
  ab:   f3 c3                   repz ret
```

```cpp-objdump:clang.txt
   0:   53                      push   rbx
   1:   48 85 d2                test   rdx,rdx
   4:   0f 84 a4 00 00 00       je     ae <_Z6NormalPK7Vector4PS_m+0xae>
   a:   4c 8d 42 ff             lea    r8,[rdx-0x1]
   e:   31 c9                   xor    ecx,ecx
  10:   48 89 c8                mov    rax,rcx
  13:   48 c1 e0 05             shl    rax,0x5
  17:   4c 8d 0c 07             lea    r9,[rdi+rax*1]
  1b:   4c 8d 54 07 10          lea    r10,[rdi+rax*1+0x10]
  20:   4c 8d 1c 06             lea    r11,[rsi+rax*1]
  24:   48 8d 44 06 10          lea    rax,[rsi+rax*1+0x10]
  29:   31 db                   xor    ebx,ebx
  2b:   0f 1f 44 00 00          nop    DWORD PTR [rax+rax*1+0x0]
  30:   39 d9                   cmp    ecx,ebx
  32:   74 65                   je     99 <_Z6NormalPK7Vector4PS_m+0x99>
  34:   66 41 0f 10 01          movupd xmm0,XMMWORD PTR [r9]
  39:   66 41 0f 10 0a          movupd xmm1,XMMWORD PTR [r10]
  3e:   66 0f 28 d0             movapd xmm2,xmm0
  42:   f2 0f 59 d2             mulsd  xmm2,xmm2
  46:   66 0f 28 d8             movapd xmm3,xmm0
  4a:   66 0f c6 db 01          shufpd xmm3,xmm3,0x1
  4f:   f2 0f 59 db             mulsd  xmm3,xmm3
  53:   f2 0f 58 da             addsd  xmm3,xmm2
  57:   66 0f 28 d1             movapd xmm2,xmm1
  5b:   f2 0f 59 d2             mulsd  xmm2,xmm2
  5f:   66 0f 28 e1             movapd xmm4,xmm1
  63:   66 0f c6 e4 01          shufpd xmm4,xmm4,0x1
  68:   f2 0f 59 e4             mulsd  xmm4,xmm4
  6c:   f2 0f 58 e2             addsd  xmm4,xmm2
  70:   f2 0f 58 e3             addsd  xmm4,xmm3
  74:   66 41 0f 10 13          movupd xmm2,XMMWORD PTR [r11]
  79:   0f 16 e4                movlhps xmm4,xmm4
  7c:   66 0f 5e c4             divpd  xmm0,xmm4
  80:   66 0f 58 c2             addpd  xmm0,xmm2
  84:   66 0f 10 10             movupd xmm2,XMMWORD PTR [rax]
  88:   66 0f 5e cc             divpd  xmm1,xmm4
  8c:   66 0f 58 ca             addpd  xmm1,xmm2
  90:   66 41 0f 11 03          movupd XMMWORD PTR [r11],xmm0
  95:   66 0f 11 08             movupd XMMWORD PTR [rax],xmm1
  99:   48 ff c3                inc    rbx
  9c:   48 39 da                cmp    rdx,rbx
  9f:   75 8f                   jne    30 <_Z6NormalPK7Vector4PS_m+0x30>
  a1:   4c 39 c1                cmp    rcx,r8
  a4:   48 8d 49 01             lea    rcx,[rcx+0x1]
  a8:   0f 85 62 ff ff ff       jne    10 <_Z6NormalPK7Vector4PS_m+0x10>
  ae:   5b                      pop    rbx
  af:   c3                      ret
```

となってました。

注目すべきは、真ん中あたりの``je     99 <_Z6NormalPK7Vector4PS_m+0x99>``から次の``jne    30 <_Z6NormalPK7Vector4PS_m+0x30>``あたりのxmmレジスタを操作してる辺り。
この辺りがコード上で処理本体になる

```cpp
const double x0 = x[i].data[0];
const double x1 = x[i].data[1];
const double x2 = x[i].data[2];
const double x3 = x[i].data[3];
const double r2 = (x0*x0 + x1*x1) + (x2*x2 + x3*x3);

const double f0 = f[i].data[0] + x0/r2;
const double f1 = f[i].data[1] + x1/r2;
const double f2 = f[i].data[2] + x2/r2;
const double f3 = f[i].data[3] + x3/r2;
f[i].data[0] = f0;
f[i].data[1] = f1;
f[i].data[2] = f2;
f[i].data[3] = f3;
```

に対応していますが、**gccはほぼ*sd命令で単発実行**（``movapd``を除く）なのに対して、**clangは*pd命令も使ってSSE2でSIMD化している**ところです。
例えば``x/r2``は、gccだとdivsdを4回で忠実に書いてますが、clangはdivpdを2回にパックしてます。
よくよく見てみると、*pdになっている部分は、後段の``f``の計算部分です。ここがSIMDで2要素同時実行できているので、実行時間は半分が半分になって3/4（つまり1.5倍の速度）という感じになっているようです。

ほんとは前段の``r2``の計算も自動ベクトル化（SIMD化）して``mulpd``がほしかったんですが、うまくコードが出てきませんでした。うーん・・・。

まぁとにかく、**clangを使うと自動ベクトル化できて、gccだとできない**という例でした。
先日の件も考えると、特にこだわりがないならgccよりclang使ったほうが良さそうですね。