4. SIMD

4.1 SIMDとは

SIMD(Single Instruction Multiple Data)とは複数のデータをベクトルレジスタに格納し一度の命令で複数の演算を行うものです。
表4-1のように現在のx86系統のCPUは標準でこの機能を持っています。
SSEでは4個、AVXでは8個の単精度の浮動小数点演算を一度に行うことができます。 従って理想的な条件ではそれぞれ4倍/8倍速くなります。 (倍精度はそれぞれ2倍/4倍)

表4-1 SIMDの分類
SIMDの種類CPUのコードネームビット数備考リリース年
SSENehalem以前128SSEからSSE4.2まで
AVXSandyBridge以降256(一部128ビット)浮動小数点演算の256ビット化2011年
AVX2Haswell以降256整数演算の256ビット化2013年

SIMD演算の代表例を図4-1に示します。
128ビットのベクトルレジスタに4個の単精度浮動小数点が入ります。 二つのベクトルレジスタ間の演算を一度に行い、結果がベクトルレジスタに入ります。


図4-1 SIMD演算例(パックド単精度浮動小数点の操作、[1]上巻図10-5)

4.2 SIMDプログラミング

SIMD演算のプログラミングには、以下の二通りがあります。
(1)アセンブラを使用する。
(2)組み込み関数を使用する。
ここではコーディングのしやすさから(2)の方法を用います。
アセンブラの命令コードと組み込み関数がほぼ1対1に対応します。 両者の計算速度は同じです。

ヘッダーファイル
組み込み関数を使用するには関数が宣言されたヘッダーファイルをincludeする必要があります。
ヘッダーファイルはSIMDの初代のMMX(64ビット演算)の"mmintrin.h"から世代を重ね、現在は"immintrin.h"です。
ヘッダーファイルは旧世代の関数を含みますので、"#include <immintrin.h>" と宣言するだけですべてのSSE/AVX関数を使用することができます。

組み込み関数
よく使う組み込み関数を表4-2に示します。 関数名はすべて"_mm"で始まります。AVXでは"_mm"が"_mm256"に変わります。
関数名の"_ps"の"p"は"Packed"を意味します。 これはレジスタの4個すべての演算を行うことを意味します。 最下位のみの演算を行うときは"s"となります。
"_ps"の"s"は単精度(float)を意味します。倍精度(double)のときは"d"となります。 倍精度では128ビットレジスタに2個のデータが入ります(AVXでは4個)。
レジスタの型はSSEでは"__m128"、AVXでは"__m256"です。 "_"は2個、"m"は1個であることに注意して下さい。
これ以外の関数については[1][2]を参考にして下さい。

表4-2 組み込み関数の一部
関数名機能対応する命令コードとその意味
_mm_load_psメモリーからレジスタに移動するMOVAPS
Move Aligned Packed Single-Precision Floating-Point Values
_mm_store_psレジスタからメモリーに移動する同上
_mm_add_psレジスタ同士の加算ADDPS
Add Packed Single-Precision Floating-Point Values
_mm_sub_psレジスタ同士の減算SUBPS
Subtract Packed Single-Precision Floating-Point Values
_mm_mul_psレジスタ同士の乗算MULPS
Multiply Packed Single-Precision Floating-Point Values

4.3 SIMDプログラミング例

リスト4-1にベクトルの内積をSIMDを用いて計算するプログラムを示します。


リスト4-1 SIMDプログラム(sdot_simd.c)
     1	/*
     2	dot product (SIMD)
     3	
     4	VC++ : cl /Ox /wd4752 sdot_simd.c
     5	gcc  : gcc -O3 -std=c99 -mavx sdot_simd.c -o sdot_simd
     6	
     7	Usage:
     8	> sdot_simd <num> <loop> [sse|avx]
     9	*/
    10	
    11	#include <stdlib.h>
    12	#include <stdio.h>
    13	#include <string.h>
    14	#include <time.h>
    15	#include <immintrin.h>
    16	
    17	static float sdot(int simd, int n, const float *a, const float *b)
    18	{
    19		int    ne = 0;
    20		float  sum = 0;
    21	
    22		if      (simd == 1) {
    23			// SSE
    24	#ifdef _WIN32
    25			__declspec(align(16)) float fsum[4];
    26	#else
    27			__attribute__((aligned(16))) float fsum[4];
    28	#endif
    29			__m128 vsum;
    30			ne = (n / 4) * 4;
    31			vsum = _mm_setzero_ps();
    32			for (int i = 0; i < ne; i += 4) {
    33				vsum = _mm_add_ps(vsum,
    34				       _mm_mul_ps(_mm_load_ps(a + i),
    35				                  _mm_load_ps(b + i)));
    36			}
    37			_mm_store_ps(fsum, vsum);
    38			sum = fsum[0] + fsum[1] + fsum[2] + fsum[3];
    39		}
    40		else if (simd == 2) {
    41			// AVX
    42	#ifdef _WIN32
    43			__declspec(align(32)) float fsum[8];
    44	#else
    45			__attribute__((aligned(32))) float fsum[8];
    46	#endif
    47			__m256 vsum;
    48			ne = (n / 8) * 8;
    49			vsum = _mm256_setzero_ps();
    50			for (int i = 0; i < ne; i += 8) {
    51				vsum = _mm256_add_ps(vsum,
    52				       _mm256_mul_ps(_mm256_load_ps(a + i),
    53				                     _mm256_load_ps(b + i)));
    54			}
    55			_mm256_store_ps(fsum, vsum);
    56			sum = fsum[0] + fsum[1] + fsum[2] + fsum[3]
    57			    + fsum[4] + fsum[5] + fsum[6] + fsum[7];
    58		}
    59	
    60		for (int i = ne; i < n; i++) {
    61			sum += a[i] * b[i];
    62		}
    63	
    64		return sum;
    65	}
    66	
    67	int main(int argc, char **argv)
    68	{
    69		int    simd = 0;
    70		int    n = 1000;
    71		int    nloop = 1000;
    72		size_t size;
    73		float  *a, *b;
    74	
    75		// arguments
    76		if (argc >= 3) {
    77			n = atoi(argv[1]);
    78			nloop = atoi(argv[2]);
    79		}
    80		if (argc >= 4) {
    81			if      (!strcmp(argv[3], "sse")) {
    82				simd = 1;
    83			}
    84			else if (!strcmp(argv[3], "avx")) {
    85				simd = 2;
    86			}
    87		}
    88	
    89		// alloc
    90		if      (simd == 1) {
    91			size = ((n + 3) / 4) * 4 * sizeof(float);
    92			a = (float *)_mm_malloc(size, 16);
    93			b = (float *)_mm_malloc(size, 16);
    94		}
    95		else if (simd == 2) {
    96			size = ((n + 7) / 8) * 8 * sizeof(float);
    97			a = (float *)_mm_malloc(size, 32);
    98			b = (float *)_mm_malloc(size, 32);
    99		}
   100		else {
   101			size = n * sizeof(float);
   102			a = (float *)malloc(size);
   103			b = (float *)malloc(size);
   104		}
   105	
   106		// setup problem
   107		for (int i = 0; i < n; i++) {
   108			a[i] = i + 1.0f;
   109			b[i] = i + 1.0f;
   110		}
   111	
   112		// timer
   113		clock_t t0 = clock();
   114	
   115		// calculation
   116		double sum = 0;
   117		for (int loop = 0; loop < nloop; loop++) {
   118			sum += sdot(simd, n, a, b);
   119		}
   120	
   121		// timer
   122		clock_t t1 = clock();
   123		double cpu = (double)(t1 - t0) / CLOCKS_PER_SEC;
   124	
   125		// output
   126		double exact = (double)nloop * n * (n + 1) * (2 * n + 1) / 6.0;
   127		printf("n=%d, nloop=%d, dot=%.6e(%.6e), cpu[sec]=%.3f\n",
   128			n, nloop, sum, exact, cpu);
   129	
   130		// free
   131		if (simd) {
   132			_mm_free(a);
   133			_mm_free(b);
   134		}
   135		else {
   136			free(a);
   137			free(b);
   138		}
   139	
   140		return 0;
   141	}

アラインメント
load/storeについては34,35,37行目または52,53,55行目のようにアラインメントされているときは "load"/"store"となりますが、アラインメントされていないときは"loadu"/"storeu"となります。
ここでアラインメントとはアドレスが16バイト(AVXでは32バイト)の倍数になるように、 コンパイラに指示することを意味します。
配列をアラインメントするには91-93行目または96-98行目のように"malloc"の代わりに"_mm_malloc"を使用し、 132,133行目のように"free"の代わりに"_mm_free"を使用します。
局所変数をアラインメントするにはVC++では25,43行目のように宣言し、 gccでは27,45行目のように宣言します。
アラインメントはメモリーアクセスを高速化するためのものですが、 アラインメントにより必ずしも速くなるとは限りません。 本ケースでも計算時間に差はないようです。
アラインメントするとコードがやや煩雑になりますので、必要に応じて使うかどうか決めて下さい。

ループ長が4または8の倍数でないとき
通常のアプリケーションではループ長に制限はなく、4または8の倍数とは限りません。
このようなときは、30,48行目のように4または8の倍数までをSIMDで計算し、 残りを60-62行目のようにSIMDを使わずに計算します。
SIMDがないときは"ne=0"となっていますので60-62行目はループ全体を表します。
以上の方法を用いるとSIMDコードと非SIMDコードを共存させることができます。
なお、20行目では和をとる変数をfloatとしていますが精度上はdoubleが望ましいです。 計算時間もほとんど変わりません。

コンパイル・リンク方法
VC++ではコンパイル・リンクオプションは特に必要ありませんが、 warningを消すには"/wd4752"を付けて下さい。 warning通り"/arch:AVX"を付けると実行プログラムはAVXをサポートしないCPUでは動きません。
> cl /Ox /wd4752 sdot_simd.c
gccではSSEを使うにはコンパイルオプションは必要ありませんが、 AVXを使うにはコンパイルオプションに"-mavx"が必要です。
> gcc -O3 -std=c99 -mavx sdot_simd.c -o sdot_simd

プログラムの実行方法
プログラムの実行方法は以下の通りです。
> sdot_simd 配列の大きさ 繰り返し回数 [sse|avx]
例えば以下のようになります。
> sdot_simd 1000000 1000 sse
最後の引数(sseまたはavx)がないときはSIMDを使用しません。
繰り返し回数は計算時間の測定誤差を小さくするためです。

4.4 SIMDの計算時間

表4-3に計算時間を示します。SIMDを使わないときとSSEまたはAVXを使ったときの結果です。
配列の大きさ(=N)と繰り返し回数(=L)の積は一定(=10^10)です。従って全体の演算量は同じです。
表の下になるほど使用メモリーが小さくなりキャッシュに乗りやすくなります。
SIMDを使わないときはNo.1-4で計算時間はほぼ一定です。
SIMDを使うときは、Windows環境ではSSEが速くなり、Linux環境ではAVXが速くなります。 他のCPUでは違う結果も得られており原因は不明です。

表4-3 SIMDの計算時間(()内はSIMDなしとの速度比)
No.配列の大きさN繰り返し回数LWindowsLinux
SIMDなしSSEAVXSIMDなしSSEAVX
110,000,0001,0008.70秒(1.0)6.37秒(1.37)6.34秒(1.37)10.46秒(1.0)7.46秒(1.40)7.06秒(1.48)
21,000,00010,0007.92秒(1.0)2.38秒(3.33)3.72秒(2.13)9.99秒(1.0)7.18秒(1.39)4.26秒(2.35)
3100,000100,0007.72秒(1.0)1.94秒(3.98)3.57秒(2.16)8.88秒(1.0)6.70秒(1.33)2.58秒(3.45)
410,0001,000,0007.72秒(1.0)1.93秒(4.00)3.56秒(2.17)9.00秒(1.0)6.72秒(1.34)1.99秒(4.52)