ベクトルの内積の計算を AVX-512 とマルチスレッドで爆速にする件。

爆速は煽りです（汗；すみません。

# 概要

AVX-512 と 16 vCPU が使える環境が手元にあったので、これらを使うとベクトルの内積をどれぐらい早く計算できるようになるのか調査してみました。

# 環境
AWS の c5.4xlarge インスタンス。8コアをハイパースレッディング（HT）で16CPUに見せるやつです。

```
$ lscpu
アーキテクチャ:                      x86_64
CPU 操作モード:                      32-bit, 64-bit
バイト順序:                          Little Endian
CPU:                                 16
オンラインになっている CPU のリスト: 0-15
コアあたりのスレッド数:              2
ソケットあたりのコア数:              8
ソケット数:                          1
NUMA ノード数:                       1
ベンダー ID:                         GenuineIntel
CPU ファミリー:                      6
モデル:                              85
モデル名:                            Intel(R) Xeon(R) Platinum 8124M CPU @ 3.00GHz
ステッピング:                        3
CPU MHz:                             3426.644
BogoMIPS:                            6000.00
ハイパーバイザのベンダー:            KVM
仮想化タイプ:                        完全仮想化
L1d キャッシュ:                      32K
L1i キャッシュ:                      32K
L2 キャッシュ:                       1024K
L3 キャッシュ:                       25344K
NUMA ノード 0 CPU:                   0-15
フラグ:                              fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ss ht syscall nx pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid aperfmperf tsc_known_freq pni pclmulqdq ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand hypervisor lahf_lm abm 3dnowprefetch invpcid_single pti fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx avx512f avx512dq rdseed adx smap clflushopt clwb avx512cd avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves ida arat pku ospke
```

# 結果

1G(1073741824)要素のベクトル同士の内積を計算するのにかかった時間。

表の読み方は以下のとおりです。

* 単位はマイクロ秒です。
* S はシングルスレッドを表します。
* SX はシングルスレッドで AVX-512 を使うことを表します。
* M はマルチスレッドを表します。
* MX はマルチスレッドで AVX-512 を使うことを表します。


## 要素が倍精度の時

```

S  	: 1400687
SX 	: 713885

スレッド数: 2
M  	: 715902
MX 	: 373578

スレッド数: 4
M  	: 372259
MX 	: 192190

スレッド数: 6
M  	: 254607
MX 	: 135453

スレッド数: 8
M  	: 191803
MX 	: 108880

スレッド数: 10
M  	: 166162
MX 	: 138235

スレッド数: 12
M  	: 139533
MX 	: 129550

スレッド数: 14
M  	: 120472
MX 	: 115902
```

## 要素が単精度の時

```
S  	: 1365853
SX 	: 356611

スレッド数: 2
M  	: 681970
MX 	: 187103

スレッド数: 4
M  	: 341022
MX 	: 96434

スレッド数: 6
M  	: 226984
MX 	: 67328

スレッド数: 8
M  	: 169235
MX 	: 54695

スレッド数: 10
M  	: 135869
MX 	: 74842

スレッド数: 12
M  	: 113181
MX 	: 67617

スレッド数: 14
M  	: 97248
MX 	: 57312
```
## 補足

* AVX-512 で最速になるように、データは 64 バイトアラインしてあり、かつデータのバイト数は64の倍数です。

* プログラムはオプション `-std=c++17 -mavx512f -Os -lpthread` でコンパイルしました。



## 感想と予測

* この章、もしかしたら間違えたこと言うかもしれないです。
* 実験は vCPU 数（16）未満のスレッド数で行ないました。

### 全体

AVX-512 とマルチスレッドを組み合わせることで

* 倍精度で 13倍弱
* 単精度で 25倍弱

の時間効率を得ました。両方とも

* 実コア数（8）スレッドのあたり

でです。

### マルチスレッドに関して

単精度はスレッド数に応じてリニアに時間効率があがります。
不思議なことに倍精度は単精度に比べて少し時間効率／スレッド数が悪いようです。
予想としては、どちらもリニアに時間効率があがる、と思っていたのですが・・原因不明です。

### AVX-512に関して

AVX-512は倍精度で8個、単精度で16個の浮動小数点演算を同時に行うことができますが、実際にはシングルスレッドで実行した場合の時間効率は以下のとおりでした。

* 倍精度で2倍弱
* 単精度で4倍弱


思ったより効率が悪いのですが、考えられる理由を下記します。

* レジスタへのローディングなどのオーバーヘッドがあるため。
* 私が何か思い切り勘違いしてる


スレッド数が8までは上記の時間効率を維持し、それ以上は時間効率がよくなりません。推測ですが、AVX-512 用のレジスタ群はコアにつき1組しかないので、コア数以上のスレッドを使ってもそれによるスピードアップは計れないのだと思われます。

Edit:
おのおののループの中で2回ずつ計算させるようにすると、それぞれ2倍近く効率アップしました。レジスタへのローディングのオーバーヘッドで間違いないようです。

# プログラム

## CPU 用

```
template	< typename T >	T
Dot( T* l, T* r, size_t size ) {
	T v = 0;
	for ( auto i = 0; i < size; i++ ) v += l[ i ] * r[ i ];
	return v;
}
```

## AVX-512 用

### 倍精度

```
double
MulAdd( double* l, double* r, size_t nBlocks, __m512d v = _mm512_setzero_pd() ) {
	for ( auto i = 0; i < nBlocks; i++ ) {
		v = _mm512_fmadd_pd( 
			_mm512_load_pd( l + i * 8 )
		,	_mm512_load_pd( r + i * 8 )
		,	v
		);
	}
	return _mm512_reduce_add_pd( v );
}
```

### 単精度

```
float
MulAdd( float* l, float* r, size_t nBlocks, __m512 v = _mm512_setzero_ps() ) {
	for ( auto i = 0; i < nBlocks; i++ ) {
		v = _mm512_fmadd_ps( 
			_mm512_load_ps( l + i * 16 )
		,	_mm512_load_ps( r + i * 16 )
		,	v
		);
	}
	return _mm512_reduce_add_ps( v );
}
```

## マルチスレッド

### CPU 用

```
template	< typename T >	T
Dot_Multi( T* l, T* r, size_t size ) {
	future< T >	wVs[ NUM_THREADS ];

	auto	wR = size % NUM_THREADS;
	auto	wD = size / NUM_THREADS;
	size_t	wIndex = 0;
	for ( auto i = 0; i < wR; i++ ) {
		wVs[ i ] = async( [=]{ return Dot( l + wIndex, r + wIndex, wD + 1 ); } );
		wIndex += wD + 1;
	}
	for ( auto i = wR; i < NUM_THREADS; i++ ) {
		wVs[ i ] = async( [=]{ return Dot( l + wIndex, r + wIndex, wD ); } );
		wIndex += wD;
	}
	T	v = 0;
	for ( auto i = 0; i < NUM_THREADS; i++ ) v += wVs[ i ].get();
	return v;
}
```

### AVX-512 用

```
template	< typename T >	T
Dot_MultiAVX512( T* l, T* r, size_t size ) {
	future< T >	wVs[ NUM_THREADS ];

	auto	nBlocks = size * sizeof( T ) / 64;

	auto	wR = nBlocks % NUM_THREADS;
	auto	wD = nBlocks / NUM_THREADS;
	size_t	wIndex = 0;
	for ( auto i = 0; i < wR; i++ ) {
		wVs[ i ] = async( [=]{ return MulAdd( l + wIndex, r + wIndex, wD + 1 ); } );
		wIndex += ( wD + 1 ) * ( 64 / sizeof( T ) );
	}
	for ( auto i = wR; i < NUM_THREADS; i++ ) {
		wVs[ i ] = async( [=]{ return MulAdd( l + wIndex, r + wIndex, wD ); } );
		wIndex += wD * ( 64 / sizeof( T ) );
	}

	T	v = 0;
	for ( auto i = nBlocks * ( 64 / sizeof( T ) ); i < size; i++ ) v += l[ i ] * r[ i ];
	for ( auto i = 0; i < NUM_THREADS; i++ ) v += wVs[ i ].get();
	return v;
}
```

# 最後に

もし上記の

* 私が何か思い切り勘違いしてる

に思い当たる点があったら教えてください。