Intel プロセッサで次の関数を実行するために、L1 キャッシュで全帯域幅を取得しようとしています。
float triad(float *x, float *y, float *z, const int n) {
float k = 3.14159f;
for(int i=0; i<n; i++) {
z[i] = x[i] + k*y[i];
}
}
これはSTREAMの triad 関数です。
この機能を備えた SandyBridge/IvyBridge プロセッサ (NASM によるアセンブリを使用) で、ピークの約 95% を取得します。ただし、Haswell を使用すると、ループを展開しない限り、ピークの 62% しか達成できません。16 回アンロールすると、92% になります。私はこれを理解していません。
NASM を使用してアセンブリで関数を記述することにしました。アセンブリのメイン ループは次のようになります。
.L2:
vmovaps ymm1, [rdi+rax]
vfmadd231ps ymm1, ymm2, [rsi+rax]
vmovaps [rdx+rax], ymm1
add rax, 32
jne .L2
Agner Fog の Optimizing Assembly マニュアルの例 12.7-12.11 で、彼はy[i] = y[i] +k*x[i]Pentium M、Core 2、Sandy Bridge、FMA4、および FMA3 に対してほぼ同じこと (ただし を除く) を行っていることがわかりました。私は彼のコードを多かれ少なかれ自力で再現することができました (実際、彼はブロードキャスト時に FMA3 の例に小さなバグを持っています)。彼は、FMA4 と FMA3 を除く各プロセッサの命令サイズ数、融合演算、実行ポートを表に示します。FMA3 用にこのテーブルを自分で作成しようとしました。
ports
size μops-fused 0 1 2 3 4 5 6 7
vmovaps 5 1 ½ ½
vfmadd231ps 6 1 ½ ½ ½ ½
vmovaps 5 1 1 1
add 4 ½ ½
jne 2 ½ ½
--------------------------------------------------------------
total 22 4 ½ ½ 1 1 1 0 1 1
サイズは命令の長さをバイト単位で示します。addand命令が μop の半分を持っている理由はjne、それらが 1 つのマクロ op に融合され (まだ複数のポートを使用する μop 融合と混同しないでください)、ポート 6 と 1 つの μop しか必要としないためです。 命令はポート 0 またはポート 1 を使用できます。ポート 0 を選択しました。負荷はポート 2 または 3 を使用できます。ポート 2 を選択し、ポート 3 を使用しました。Agner Fog のテーブルと一貫性を保つために、異なるポートに行くことができる命令は 1/2 の確率で各ポートに等しく行くと言う方が理にかなっていると思うので、ポートに 1/2 を割り当てて、vfmadd231psvmovapsvfmadd231psvmovaps行くvmadd231psことができます。に。
この表と、すべての Core2 プロセッサがクロック サイクルごとに 4 つの μops を実行できるという事実に基づくと、このループはクロック サイクルごとに可能であるように見えますが、私はそれを得ることができませんでした。アンロールしないと、Haswell でこの関数のピーク帯域幅に近づけない理由を誰か説明してもらえますか? これは展開せずに可能ですか？この関数の ILP を最大化しようとしていることをはっきりさせておきます (最大帯域幅だけが必要なわけではありません)。それが、展開したくない理由です。
編集: Iwillnotexist Idonotexist が IACA を使用して、ストアがポート 7 を使用しないことを示したので、ここに更新があります。アンロールせずに 66% バリアを破り、アンロールせずに反復ごとに 1 クロック サイクルでこれを行うことができました (理論的には)。まず、店舗の問題を解決しましょう。
[base + offset]Stephen Canon はコメントで、ポート 7 のアドレス生成ユニット (AGU) はや notなどの単純な操作しか処理できないと述べ[base + index]ました。Intel 最適化リファレンス マニュアルで見つけたのは、port7 に関するコメントで、「Simple_AGU」と書かれており、simple の意味が定義されていませんでした。しかし、その後 Iwillnotexist Idonotexist は、IACA のコメントで、この問題は 6 か月前に Intel の従業員が 2014 年 3 月 11 日に書いた記事で既に言及されていることを発見しました。
Port7 AGU は、単純なメモリ アドレス (インデックス レジスタなし) を持つストアでのみ機能します。
Stephen Canon は、「ストア アドレスをロード オペランドのオフセットとして使用する」ことを提案しています。私はこのようにこれを試しました
vmovaps ymm1, [rdi + r9 + 32*i]
vfmadd231ps ymm1, ymm2, [rsi + r9 + 32*i]
vmovaps [r9 + 32*i], ymm1
add r9, 32*unroll
cmp r9, rcx
jne .L2
これにより、実際にストアはポート 7 を使用します。ただし、vmadd231psIACA からわかるように負荷と融合しないという別の問題があります。また、cmp私の本来の機能にはなかった命令も追加で必要です。そのため、ストアで使用する micro-op は 1 つ少なくなりますが、(cmpマクロが と融合するため) にはもう 1 つ必要です。IACA は、1.5 のブロック スループットを報告しています。実際には、これはピークの約 57% しか得られません。addcmpjne
しかし、vmadd231ps命令を負荷と融合させる方法も見つけました。これは、このようにアドレス指定 [絶対 32 ビット アドレス + インデックス] を使用する静的配列を使用してのみ行うことができます。Evgeny Kluev オリジナルはこれを提案しました。
vmovaps ymm1, [src1_end + rax]
vfmadd231ps ymm1, ymm2, [src2_end + rax]
vmovaps [dst_end + rax], ymm1
add rax, 32
jl .L2
src1_end、src2_end、およびは、静的dst_end配列の終了アドレスです。
これは、私が予想した 4 つの融合マイクロオペレーションを使用して、質問の表を再現しています。これを IACA に入れると、1.0 のブロック スループットが報告されます。理論的には、これは SSE および AVX バージョンと同様に機能するはずです。実際には、ピークの約 72% を取得します。これで 66% の壁は破れましたが、16 回アンロールした 92% にはほど遠い状態です。したがって、Haswell では、ピークに近づくための唯一のオプションは展開することです。これは Ivy Bridge 経由の Core2 では必要ありませんが、Haswell では必要です。
編集終了:
これをテストするための C/C++ Linux コードを次に示します。NASM コードは、C/C++ コードの後に掲載されています。変更する必要があるのは、周波数番号だけです。行のdouble frequency = 1.3;1.3 を、プロセッサの動作 (公称ではない) 周波数に置き換えます (BIOS でターボが無効になっている i5-4250U の場合は 1.3 GHz です)。
でコンパイル
nasm -f elf64 triad_sse_asm.asm
nasm -f elf64 triad_avx_asm.asm
nasm -f elf64 triad_fma_asm.asm
g++ -m64 -lrt -O3 -mfma tests.cpp triad_fma_asm.o -o tests_fma
g++ -m64 -lrt -O3 -mavx tests.cpp triad_avx_asm.o -o tests_avx
g++ -m64 -lrt -O3 -msse2 tests.cpp triad_sse_asm.o -o tests_sse
C/C++ コード
#include <x86intrin.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#define TIMER_TYPE CLOCK_REALTIME
extern "C" float triad_sse_asm_repeat(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_sse_asm_repeat_unroll16(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_avx_asm_repeat(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_avx_asm_repeat_unroll16(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_fma_asm_repeat(float *x, float *y, float *z, const int n, int repeat);
extern "C" float triad_fma_asm_repeat_unroll16(float *x, float *y, float *z, const int n, int repeat);
#if (defined(__FMA__))
float triad_fma_repeat(float *x, float *y, float *z, const int n, int repeat) {
float k = 3.14159f;
int r;
for(r=0; r<repeat; r++) {
int i;
__m256 k4 = _mm256_set1_ps(k);
for(i=0; i<n; i+=8) {
_mm256_store_ps(&z[i], _mm256_fmadd_ps(k4, _mm256_load_ps(&y[i]), _mm256_load_ps(&x[i])));
}
}
}
#elif (defined(__AVX__))
float triad_avx_repeat(float *x, float *y, float *z, const int n, int repeat) {
float k = 3.14159f;
int r;
for(r=0; r<repeat; r++) {
int i;
__m256 k4 = _mm256_set1_ps(k);
for(i=0; i<n; i+=8) {
_mm256_store_ps(&z[i], _mm256_add_ps(_mm256_load_ps(&x[i]), _mm256_mul_ps(k4, _mm256_load_ps(&y[i]))));
}
}
}
#else
float triad_sse_repeat(float *x, float *y, float *z, const int n, int repeat) {
float k = 3.14159f;
int r;
for(r=0; r<repeat; r++) {
int i;
__m128 k4 = _mm_set1_ps(k);
for(i=0; i<n; i+=4) {
_mm_store_ps(&z[i], _mm_add_ps(_mm_load_ps(&x[i]), _mm_mul_ps(k4, _mm_load_ps(&y[i]))));
}
}
}
#endif
double time_diff(timespec start, timespec end)
{
timespec temp;
if ((end.tv_nsec-start.tv_nsec)<0) {
temp.tv_sec = end.tv_sec-start.tv_sec-1;
temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec;
} else {
temp.tv_sec = end.tv_sec-start.tv_sec;
temp.tv_nsec = end.tv_nsec-start.tv_nsec;
}
return (double)temp.tv_sec + (double)temp.tv_nsec*1E-9;
}
int main () {
int bytes_per_cycle = 0;
double frequency = 1.3; //Haswell
//double frequency = 3.6; //IB
//double frequency = 2.66; //Core2
#if (defined(__FMA__))
bytes_per_cycle = 96;
#elif (defined(__AVX__))
bytes_per_cycle = 48;
#else
bytes_per_cycle = 24;
#endif
double peak = frequency*bytes_per_cycle;
const int n =2048;
float* z2 = (float*)_mm_malloc(sizeof(float)*n, 64);
char *mem = (char*)_mm_malloc(1<<18,4096);
char *a = mem;
char *b = a+n*sizeof(float);
char *c = b+n*sizeof(float);
float *x = (float*)a;
float *y = (float*)b;
float *z = (float*)c;
for(int i=0; i<n; i++) {
x[i] = 1.0f*i;
y[i] = 1.0f*i;
z[i] = 0;
}
int repeat = 1000000;
timespec time1, time2;
#if (defined(__FMA__))
triad_fma_repeat(x,y,z2,n,repeat);
#elif (defined(__AVX__))
triad_avx_repeat(x,y,z2,n,repeat);
#else
triad_sse_repeat(x,y,z2,n,repeat);
#endif
while(1) {
double dtime, rate;
clock_gettime(TIMER_TYPE, &time1);
#if (defined(__FMA__))
triad_fma_asm_repeat(x,y,z,n,repeat);
#elif (defined(__AVX__))
triad_avx_asm_repeat(x,y,z,n,repeat);
#else
triad_sse_asm_repeat(x,y,z,n,repeat);
#endif
clock_gettime(TIMER_TYPE, &time2);
dtime = time_diff(time1,time2);
rate = 3.0*1E-9*sizeof(float)*n*repeat/dtime;
printf("unroll1 rate %6.2f GB/s, efficency %6.2f%%, error %d\n", rate, 100*rate/peak, memcmp(z,z2, sizeof(float)*n));
clock_gettime(TIMER_TYPE, &time1);
#if (defined(__FMA__))
triad_fma_repeat(x,y,z,n,repeat);
#elif (defined(__AVX__))
triad_avx_repeat(x,y,z,n,repeat);
#else
triad_sse_repeat(x,y,z,n,repeat);
#endif
clock_gettime(TIMER_TYPE, &time2);
dtime = time_diff(time1,time2);
rate = 3.0*1E-9*sizeof(float)*n*repeat/dtime;
printf("intrinsic rate %6.2f GB/s, efficency %6.2f%%, error %d\n", rate, 100*rate/peak, memcmp(z,z2, sizeof(float)*n));
clock_gettime(TIMER_TYPE, &time1);
#if (defined(__FMA__))
triad_fma_asm_repeat_unroll16(x,y,z,n,repeat);
#elif (defined(__AVX__))
triad_avx_asm_repeat_unroll16(x,y,z,n,repeat);
#else
triad_sse_asm_repeat_unroll16(x,y,z,n,repeat);
#endif
clock_gettime(TIMER_TYPE, &time2);
dtime = time_diff(time1,time2);
rate = 3.0*1E-9*sizeof(float)*n*repeat/dtime;
printf("unroll16 rate %6.2f GB/s, efficency %6.2f%%, error %d\n", rate, 100*rate/peak, memcmp(z,z2, sizeof(float)*n));
}
}
System V AMD64 ABI を使用した NASM コード。
triad_fma_asm.asm:
global triad_fma_asm_repeat
;RDI x, RSI y, RDX z, RCX n, R8 repeat
;z[i] = y[i] + 3.14159*x[i]
pi: dd 3.14159
;align 16
section .text
triad_fma_asm_repeat:
shl rcx, 2
add rdi, rcx
add rsi, rcx
add rdx, rcx
vbroadcastss ymm2, [rel pi]
;neg rcx
align 16
.L1:
mov rax, rcx
neg rax
align 16
.L2:
vmovaps ymm1, [rdi+rax]
vfmadd231ps ymm1, ymm2, [rsi+rax]
vmovaps [rdx+rax], ymm1
add rax, 32
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
global triad_fma_asm_repeat_unroll16
section .text
triad_fma_asm_repeat_unroll16:
shl rcx, 2
add rcx, rdi
vbroadcastss ymm2, [rel pi]
.L1:
xor rax, rax
mov r9, rdi
mov r10, rsi
mov r11, rdx
.L2:
%assign unroll 32
%assign i 0
%rep unroll
vmovaps ymm1, [r9 + 32*i]
vfmadd231ps ymm1, ymm2, [r10 + 32*i]
vmovaps [r11 + 32*i], ymm1
%assign i i+1
%endrep
add r9, 32*unroll
add r10, 32*unroll
add r11, 32*unroll
cmp r9, rcx
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
triad_ava_asm.asm:
global triad_avx_asm_repeat
;RDI x, RSI y, RDX z, RCX n, R8 repeat
pi: dd 3.14159
align 16
section .text
triad_avx_asm_repeat:
shl rcx, 2
add rdi, rcx
add rsi, rcx
add rdx, rcx
vbroadcastss ymm2, [rel pi]
;neg rcx
align 16
.L1:
mov rax, rcx
neg rax
align 16
.L2:
vmulps ymm1, ymm2, [rdi+rax]
vaddps ymm1, ymm1, [rsi+rax]
vmovaps [rdx+rax], ymm1
add rax, 32
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
global triad_avx_asm_repeat2
;RDI x, RSI y, RDX z, RCX n, R8 repeat
;pi: dd 3.14159
align 16
section .text
triad_avx_asm_repeat2:
shl rcx, 2
vbroadcastss ymm2, [rel pi]
align 16
.L1:
xor rax, rax
align 16
.L2:
vmulps ymm1, ymm2, [rdi+rax]
vaddps ymm1, ymm1, [rsi+rax]
vmovaps [rdx+rax], ymm1
add eax, 32
cmp eax, ecx
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
global triad_avx_asm_repeat_unroll16
align 16
section .text
triad_avx_asm_repeat_unroll16:
shl rcx, 2
add rcx, rdi
vbroadcastss ymm2, [rel pi]
align 16
.L1:
xor rax, rax
mov r9, rdi
mov r10, rsi
mov r11, rdx
align 16
.L2:
%assign unroll 16
%assign i 0
%rep unroll
vmulps ymm1, ymm2, [r9 + 32*i]
vaddps ymm1, ymm1, [r10 + 32*i]
vmovaps [r11 + 32*i], ymm1
%assign i i+1
%endrep
add r9, 32*unroll
add r10, 32*unroll
add r11, 32*unroll
cmp r9, rcx
jne .L2
sub r8d, 1
jnz .L1
vzeroupper
ret
triad_sse_asm.asm:
global triad_sse_asm_repeat
;RDI x, RSI y, RDX z, RCX n, R8 repeat
pi: dd 3.14159
;align 16
section .text
triad_sse_asm_repeat:
shl rcx, 2
add rdi, rcx
add rsi, rcx
add rdx, rcx
movss xmm2, [rel pi]
shufps xmm2, xmm2, 0
;neg rcx
align 16
.L1:
mov rax, rcx
neg rax
align 16
.L2:
movaps xmm1, [rdi+rax]
mulps xmm1, xmm2
addps xmm1, [rsi+rax]
movaps [rdx+rax], xmm1
add rax, 16
jne .L2
sub r8d, 1
jnz .L1
ret
global triad_sse_asm_repeat2
;RDI x, RSI y, RDX z, RCX n, R8 repeat
;pi: dd 3.14159
;align 16
section .text
triad_sse_asm_repeat2:
shl rcx, 2
movss xmm2, [rel pi]
shufps xmm2, xmm2, 0
align 16
.L1:
xor rax, rax
align 16
.L2:
movaps xmm1, [rdi+rax]
mulps xmm1, xmm2
addps xmm1, [rsi+rax]
movaps [rdx+rax], xmm1
add eax, 16
cmp eax, ecx
jne .L2
sub r8d, 1
jnz .L1
ret
global triad_sse_asm_repeat_unroll16
section .text
triad_sse_asm_repeat_unroll16:
shl rcx, 2
add rcx, rdi
movss xmm2, [rel pi]
shufps xmm2, xmm2, 0
.L1:
xor rax, rax
mov r9, rdi
mov r10, rsi
mov r11, rdx
.L2:
%assign unroll 8
%assign i 0
%rep unroll
movaps xmm1, [r9 + 16*i]
mulps xmm1, xmm2,
addps xmm1, [r10 + 16*i]
movaps [r11 + 16*i], xmm1
%assign i i+1
%endrep
add r9, 16*unroll
add r10, 16*unroll
add r11, 16*unroll
cmp r9, rcx
jne .L2
sub r8d, 1
jnz .L1
ret